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ABSTRACT 


In  the  first  part  of  this  thesis  solid-fluid  chemical 
reactors,  with  reaction  if  any  only  in  the  fluid  phase,  are 
simulated  mathematically  with  due  consideration  of  the  tem¬ 
perature  gradients  in  the  solid  particles.  A  heat  diffusion 
equation  describes  the  temperature  change  with  time  in  the 
solid  spheres,  a  separate  equation  describes  the  gas  tempera¬ 
ture  along  the  reactor .  They  are  coupled  together  by  the 
boundary  condition  on  the  surface  of  the  solids.  If  reaction 
is  present  additional  equations  for  component  balances  along 
the  reactor  are  also  included. 

For  a  moving  bed  reactor,  assuming  a  uniform  velocity 
of  the  bed,  the  system  can  be  readily  expressed  by  a  set  of 
simultaneous  ordinary  differential  equations,  having  z  the 
length  along  the  reactor  as  the  sole  independent  variable,  the 
dependent  variables  being  composition,  temperature  and  pressure 
of  the  fluid  and  a  series  of  temperatures  along  the  radius  of 
the  solid  particles.  It  is  then  solved  firstly  by  an  analog 
computer  for  no  chemical  reaction,  and  numerically  by  the 
Runge-Kutta-Gill  method  in  a  digital  computer  for  both  cases 
with  and  without  chemical  reactions. 

However,  in  the  fixed  bed  regenerative  reactor  the 
time  variables  for  both  the  solid  and  gas  are  different.  The 
above  approach  is  not  applicable.  A  mathematical  model  is 
proposed  in  which  the  surface  temperature  in  a  differential 


volume  of  the  reactor  is  assumed  to  change  linearly  with  time, 
so  that  the  gas  phase  in  passing  through  this  reactor  volume 
sees  the  mean  surface  temperature  of  the  solid.  Thus,  the 
equations  describing  the  gas  phase  and  the  solid  phase  are 
effectively  decoupled.  With  an  assumed  final  surface 
temperature  for  a  given  time,  an  analytical  solution  can  be 
obtained  for  the  solid  phase,*  the  gas  phase  is  described  by  a 
set  of  ordinary  differential  equations  which/with  the  mean 
surface  temperature  as  a  parameter,  can  be  readily  solved 
numerically.  The  correct  final  surface  temperature  is 
established  through  trial-and-error  using  an  energy  balance 
on  the  two  phases. 

The  above  proposed  methods  had  been  tested  numer¬ 
ically  against  known  solutions  for  simple  cases  or  against 
established  methods  for  more  general  cases  and  found  to  be 
accurate  and  feasible. 

The  latter  part  of  the  thesis  involves  the  study  of 
the  correlation  of  multivariate  data.  A  least  squares  analysis 
with  self-generating  orthogonal  polynomials  as  approximating 
functions  was  used  for  curve  fitting.  It  was  found  to  be 
most  convenient  to  use  and  gave  good  accuracy  for  high  order 
approximations.  This  approach  was  developed  to  surface  fitting, 
with  two  independent  variables.  The  results  obtained  in¬ 
dicate  the  reliability  and  the  accuracy  of  this  method. 

If  the  approximating  functions  are  of  some  specified 
form  the  polynomial  approach  will  no  longer  be  applicable. 


If  the  Chebyshev  criterion  of  fit  can  be  used,  then  the 
analysis  becomes  a  linear  programming  problem  with  the  maximum 
absolute  deviation  to  be  minimized  as  the  objective  function. 
The  problem  can  be  solved  more  readily  through  its  dual,  using 
the  Simplex  method.  A  program  for  the  Simplex  method  has  been 
developed  for  the  dual  problem  in  search  for  the  optimal  base. 
Thus,  a  corresponding  optimal  in  the  primal  can  be  obtained, 
which  is  a  set  of  equations  equal  to  the  number  of  coefficients 
plus  the  minimal  maximum  deviation.  This  may  be  solved  by 
any  conventional  method. 

It  was  applied  to  curve  fitting  and  the  evaluation 
of  Benedict-Webb-Rubin  equation  of  state.  The  results  of  the 
approximation  illustrate  the  accuracy  of  the  method. 
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PART  ONE 

COMPUTER  MODELS  OF  SOLIDS-FLUID  CHEMICAL  REACTORS 


REACTIONS  IN  FLUID  PHASE  ONLY 
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1 •  INTRODUCTION 

Solids-gas  contacting  operations  are  among  the  most 
important  and  widely  encountered  in  the  chemical  metallurgical 
and  petroleum  industries.  A  detailed  description  of  this 
operation  is  given  by  Vener  and  Robinson(l) .  From  the  point 
of  view  of  chemical  engineering,  solids-gas  contacting  should 
be  considered  as  a  broad  unit  operation,  in  which  moving-bed, 
fixed-bed  and  f luidized-bed  techniques  may  be  considered  as 
some  typical  examples. 

In  the  processes  of  adsorption,  drying  or  solvent 
removal,  thermal  decomposition  of  solids,  thermal  cracking 
and  pebble-bed  heat  exchangers,  the  temperature  or  concentra¬ 
tion  gradients  in  the  solids  are  of  prime  concern  in  the  design 
of  the  reactors.  It  is  the  purpose  of  this  thesis  to  study 
the  temperature  or  concentration  profiles  along  the  reactors 
at  various  times  of  the  operation  with  due  consideration  of  the 
temperature  or  concentration  gradients  in  the  solids.  How¬ 
ever,  in  the  fluidized  bed  the  solid  diameters  are  so  small 
that  the  temperature  gradients  in  the  solids  can  probably  be 
neglected.  Thus,  only  moving  beds  and  fixed  bedswill  be 
studied  in  the  following  sections. 

A  moving  bed  can  be  defined  as  a  body  of  solid 
particles  moving  downward  by  gravity  in  a  reactor,  through 
which  gas  may  flow  parallel  or  countercurrent  to  the  solids. 

Due  to  its  ability  to  handle  a  wide  variety  of  solids, 
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elimination  of  dust  entrainment  problems,  reasonable  pressure 
drops,  simple  provisions  for  establishing  desirable  temperature 
gradients  along  the  reactor  to  permit  optimum  conversions,  uni¬ 
formity  in  residence  time  of  solids,  suitable  temperature  con¬ 
trol  by  various  heat  transfer  techniques  and  many  other  advan¬ 
tages,  the  moving  bed  process  is  considered  to  have  one  of  the 
best  growth  potentials  in  the  solids-gas  contacting  field(2). 

The  process  of  bringing  gases  into  contact  with 
fixed  solids  is  even  more  common.  A  fixed  bed  may  be  defined 
as  a  body  of  solid  particles  packed  in  a  reactor  through  which 
gases  may  pass  upward  or  downward.  As  compared  to  the  moving 
bed,  the  main  advantage  of  fixed  bed  processes  is  the  lack 
of  necessity  of  equipment  for  circulating  the  solids.  How¬ 
ever,  owing  to  the  relatively  poor  heat  transfer  character¬ 
istics,  excessive  temperature  gradients  often  exist  within 
the  bed . 

The  problems  to  be  considered  involve  the  calcula¬ 
tion  of  the  behavior  of  a  flowing  gas  stream  and  a  moving 
bed  or  fixed  bed  of  solids,  where  chemical  reaction, if  any 
occurs  only  in  the  gas  phase.  The  analysis  is  predicated  upon 
the  following  assumptions:  no  radial  temperature  or  concentration 
gradients  in  the  gas  phase,  i.e.  one-dimensional  reactor;  no 
axial  dispersion  of  either  energy  or  mass  except  by  bulk  flow; 

V.he  solids  are  homogeneous  uniform-sized  spheres;  the  reactors 
are  adiabatic;  and  the  particles  are  large  enough  so 
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that  internal  temperature  or  concentration  gradients  are 
significant  but  small  enough  so  that  the  fluid  environment  can 
be  considered  as  uniform,  thus  ensuring  radial  symmetry. 
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II,  STUDIES  OF  FINITE  DIFFERENCE  TECHNIQUES  IN 

HEAT  TRANSFER  BETWEEN  SOLIDS  AND  GAS 


Before  attempting  to  solve  the  design  problem,  it  is 
appropriate  to  consider  the  techniques  which  are  presently 
used  in  simple  systems. 

Consider  a  sphere  of  homogeneous  and  isotropic  solid, 
having  initial  temperature  TQ.  It  is  submerged  at  time  zero 
into  a  well-stirred  fluid  of  constant  temperature,  Tg.  The 
temperature  distribution  in  the  sphere  at  any  time  is  required. 

Mathematically 


St 

St 


K  ( 


S2T 

Sr^ 


+ 


2  ST 

r  Sr 


(II. 1) 


with  boundary  conditions  of 


ST 

Sr 


=  0 


r  =  0 


St 

Sr 


hf 


(Tg  "  Ts) 


r  =  R 


(II. 2) 


T  =  T 
s 


and  initial  condition 


r  =  R 


T  =  T, 


t  =  0  for  0 <  r  <  R  (II. 3) 


This  system  is  to  be  solved  by  three  different 
methods,  namely,  1)  analytical,  2)  implicit  and  3)  explicit 


method . 
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1 )  Analytical  Method : 


By  separation  of  variables,  equation  (II.l)  can  be 
readily  solved  to  give 


oo 


T  = 


T  + 

g 


n=l  v'  ^ 


{  9  'i  A_  sin  a„r 

| exp  (-Kan  t)|  - 


(II. 4) 


where  a^  are  the  roots  of  the  transcendental  equation 

(anR)  ks 

tan  (anR)  =  - -  (II.  5) 

k  -  Rh  - 
s  f  * 


The  constants  A^  are  the  coefficients  of  the  orthogonal 
functions  in  approximating  the  initial  temperature  distribu¬ 
tion  and  are  given  as 


(sin  anr)/r 


(Tq  -  Tg)  dr 


R 


o 


anw 


( sin 


v  2  /  2 

anr)  /r 


dr 


sin  anR  -  anR  cos  a  R 


anR  -  1/2  sin  2anR  . 


(II. 6) 


2 )  Implicit  (Crank-Nicholson) Method : 

According  to  the  implicit  expression  proposed  by 
Crank-Nicholson ( 3) ,  the  finite  difference  representation  of 
equation  (II.l)  according  to  the  network  given  in  Figure  la, 
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is 


(6W  +  1)^  -  6WT2 


k+1 


(1  -  6W)T]_  -I-  6WT2 


WT 


1  +  (1  +  2W/3)T2  -  (5W^Jl*2 


J  k+1 


-WT 


x  +  (1  -  2W/3)T2  +(5W/3)T3 


•W j  '  /3  T  j_2  -  W(  1  -  2  j  '  )  T  j_1  +  (1  -  W(j'  -  2))Tj 


-  W(1  +  2 j 1 /3) T . 


j+1 


k+1 


Wj'/3  T._9  +  W(1  -  2 j ' ) T .  ,  +  (1  +  W ( j  —  2 ) ) T • 
J  ^  J~+  J 


+  W(1  +  2 j ‘/3)T 


j  +  1 


for  j  =  3#  4,  5,  ...  ,  M-l 


(II.  7) 


(-WM'/3  -  2c)Tm_2  -  ( W( 1  -  2M ' )  -  9C)  TM_1 


+  (1  -  W(M '  -  2)  -  18C) T 


M 


k+1 


(WM'/3  +  2C)Tm_2  +  (W(l  -  2M‘)  -  9C)  TM-1 


+  (1  +  W(M '  -  2)  +  18C)  T 


M 


k  +  2BCTg 


where  TM+1  =  (2TM_2  -  9TM-1  +  18TM  +  BT  )/(ll  +  B)  (II. 8) 


and 


At  *K 


W  = 


2  -Ar2 


B  = 


12  Ar  hf 


3  = 


(j  -  1) 


M*  =  1/(M-1), 


C  =  W(1  +  2M'/3)/(ll  +  B) 


Expression  (II. 7)  is  the  final  computational  system.  Other 
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Figure  1 


Network  for  Heat  Conduction  Equation. 

— rm»  vm-anmi  ^Tnnrt  ,  xw i  ■  »u  i  ■  i  i  ,  »i  »  1 1  up  m  , , ,  ■■■■—.  — _  .  _ _ 


(a) 


T1  t2 


TM+1 


T  j  t  k  =  T(jAr,kat) 


(b) 


« - — e - « - « - e «  Tc 

O 


T  .  ,  =  T  (  ( j-^)ar,kAt) 

J  t  ^ 


tM-1  TM 
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forms  of  implicit  expressions  are  possible  and  are  discussed 

in  detail  by  Lapidus(4) . 

3)  Explicit  (Runge-Kutta-Gill)  Method: 

The  original  form  of  this  numerical  integration  of 

ordinary  differential  equations  was  due  to  Runge(5),  but  was 

later  modified  by  Kutta(6) .  A  complete  derivation  of  the 

method  is  very  complicated  and  can  be  found  in  a  paper  by 

Ince(7).  The  method  has  been  further  refined  for  digital 

computation  by  Gill (8)  so  as  to  give  better  accuracy  and  to 

minimize  storage  requirements. 

Here  R  is  again  divided  into  M  equal  sections,  with 

discrete  functional  values  of  T  assumed  to  be  in  the  middle 

of  each  section,  representing  the  mean  value  of  the  spherical 

shell  as  sketched  in  Figure  lb.  Now  third  order  correct 

2  2 

expressions  are  used  for  dT/dr  and  o  T/dr  as  derived  in  . 
Appendix  A.  Thus,  equation  (II .1)  and  its  associated  boundary 
conditions  give  rise  to  a  set  of  simultaneous  first  order 
ordinary  differential  equations  as 


K 


DT- 


8Ar‘ 


( -26T-,  +  27T0  -  T~) 


K 


DT, 


9Ar 


-  ( -T,  -  12T0  +  13T-) 

2  1  2  3 


K 


DT 


7  (77-7  Ti-: 


2j-5 


j-2  +  2 j-1  i^_1  2 j-1  '  6 j-3  ^+1 


4j-4  6j+l 

T_i  +  -  T  •  ,  ,  ) 


Ar  6  j  -  3 


j  =  3 ,  4 ,  . .  . ,  M- 1 


* 
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6 


-  (97  -  12M  +  18) T 


M-2 


(50y  -  120M  +  140)  T 


+  (225Y  -  300M  +  90)Tm  +  0.?  -  T  |  (II. 9) 

where 

D  =  d/dt 
f3  =  60Arh^/kg 

y  =  ( 192M  +  32 ) /( 184  +  0) 

<5  =  K/ ( Ar 2 • 30 ( 2M- 1 ) ) 


However  at  r  =  R 


dT 


s 


dr 


=  (184Tg  -  225Tm  +  50Tm_1  -  9TM_2 ) /( 60Ar ) 

T  =  (£T  +  22 5T  -  50T  n  +  9T  0)/(184  +  p)  (II. 10) 
s  g  M  M-l  M-2 

System  (II. 9)  is  then  integrated  numerically  by  the  Runge- 
Kutta-Gill  method  which  is  described  inAppendix  B.  It  can 
be  seen  that  the  computational  program  for  the  explicit  formula 
(II. 9)  is  quite  simple  and  direct  in  comparison  with  that  of 
the  implicit  formula  (II. 7).  However,  there  are  some  severe 
computational  restrictions  involved  in  the  use  of  the  formula. 
The  stability  requirement  may  demand  such  small  time  increments 
that  computing  times  become  excessive.  An  analysis  of  the 
stability  is  presented  in  Appendix  C. 
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II. 1  Numerical  Example 

For  numerical  evaluation,  the  following  specific 


values  are  used: 

T 

o 

0 

R  = 

0.1  (ft) 

hf  = 

50  BTU/(hr) (ft2) 

^s  = 

2.0  BTU/(hr) (ft2) ( °F/ft) 

T  = 

g 

1.0 

Ar  = 

0.01  (ft) 

K  = 

0.01  ( ft2)  /(hr) 

The  results  from  the  above  three  methods  are  plotted  in 
Graph  No.  1,  which  shows  that  the  Runge-Kutta  method  gives 
slightly  better  accuracy  than  that  of  the  implicit  Crank - 
Nicholson  method.  This  difference  in  the  results  is  expected 
as  the  finite  difference  representation  of  dT/dt  is  fourth 
order  correct  for  the  Runge-Kutta  method,  but  only  second 
order  correct  for  the  implicit  method.  However,  the 
implicit  method  cannot  be  conveniently  used  if  any  of  the 
differential  equations  are  non-linear. 


Normalized  Temperature  of  Solid 
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Graph  No.  1 

Comparison  of  Results  by  Various  Finite  Difference 

Methods  in  Heat  Conduction. 


Analytical  solution 
Explicit  (R-K-G)  method 


Normalized  Distance  along  the  Radius 
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III«  MOVING  BED  CHEMICAL  REACTOR  DESIGN 


For  a  differential  section  of  the  reactor,  the 
energy  balance  in  the  gas  phase  is 


dT. 


-4(hw/DR>  (T 


-  T  )  ~ 
w' 


hfk(T  -  TJ 
r  s  g  s 


R(x,  P ,  Tg)  AH^, 


dz 


(III.l) 


The  fractional  conversion  of  the  limiting  reactant  in  the 
flow  reactor  is 


dx 

dz 


R  (x, P , Tg) 


(III. 2) 


and  the  pressure  drop  is 


dP 

dz 


1.635  G  (1  -  6  ) 


(III. 3) 


(  1 . 7  7  D  -  (1  -  €  )) 


0.05 


-  °-81 j  gc  Dp  €  Pg 


which  is  proposed  by  Ergun  and  modified  by  Glaser  and  Thodos(9) 
to  take  into  account  the  motion  of  the  solid  particles.  If 
there  are  more  than  one  independent  reaction  then  additional 
equations  of  the  form  of  equation  (III. 2)  must  be  included. 

Both  axial  and  radial  diffusions  have  been  neglected.  These 
then  are  the  equations  describing  a  one-dimensional  reaction 
in  which  z  is  the  independent  variable  and  Tg/  X  and  P  the 
dependent  variables.  The  only  term  in  the  equations,  which 
cannot  be  immediately  expressed  in  terms  of  Tg,  X  and  P,  is  Ts, 
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the  surface  temperature  of  the  moving  bed  of  solids.  To 
evaluate  it,  an  energy  balance  must  be  written  about  the  solid- 
gas  interface. 

Therefore,  attention  must  now  be  focussed  on  the 
solid  phase,  where  the  thermal  behavior  of  a  homogeneous  solid 
sphere  is  given  by 


St  d2T 

—  =  K  (— 


St 


Sr' 


2  St 

+ - ) 

r  Sr 


(III. 4) 


where  K,  the  thermal  diffusivity  (k  _/p_c  )  is  assumed  to  be 

O  b  D 

independent  of  r  and  T. 

In  the  problem  under  study,  since  the  solids  are 
moving  through  the  reactor  at  a  constant  rate,  t  and  z  are 
related  linearly  and  equation  (ill. 4)  may  be  written  as 


St 

Sz 


=  B  ( 


S2t 


Sr‘ 


St 


2 

+  ; 


(III. 5) 


where 


B  = 


K  p-  (1  -€  ) 


The  boundary  conditions  are 

St 


Sr 


=  0 


St  h£ 

Sr  k~ 


<Tg  '  V 


r  =  0 


r  =  R 


(III. 6) 
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and  initial  conditions  are 

=  T 

o 

=  Tg° 

at  z  =  0 

=  xo 
=  po 

The  complexity  of  the  solution  depends  upon  the 
relationship  that  exists  between  T  and  other  variables. 

o 

Furthermore,  equations  (III.1)  to  (III. 3)  are  inherently  non¬ 
linear,  In  essence,  the  problem  involves  the  solution  of  a 
second  order,  parabolic  partial  differential  equation  describ¬ 
ing  the  solid,  a  set  of  first  order  ordinary  differential 
equations  describing  the  gas,  coupled  by  a  heat  transfer 
rate  equation  (III. 6)  at  the  gas-solid  interface. 

The  proposed  solution  uses  numerical  finite  dif¬ 
ference  methods,  outlined  partly  in  section  II.  The  first 
step  is  to  convert  the  partial  differential  equation  (III  .  5) 
into  a  set  of  ordinary  equations  by  replacing  dT/dr  and 
d  T/dr  at  various  discrete  values  of  r,  with  their  appropriate 
finite-difference  expressions  from  Appendix  A.  As  recommended 
in  section  II,  the  network  in  Figure  lb  is  used  and  the  various 
equations  are  shown  in  Table  1. 

The  system  of  differential  equations  to  be  solved 
comprise  equations  (III.  1), (ill. 2)  ,  (III. 3)  and  the  set  (III. 7). 
The  numerical  solution  of  these  equations,  much  too  formidable 
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TABLE  1 


Differential  Equations  f or  Te nro o r r. t u r e 

Distribution  in  Spheres 


1  1 2 


T10 


r=0 


-Ar 


dTn  B 


dz  8 


(-26T1  +  27T2 


T3) 


dT' 


dz 


B  ' 

9 


(  -  T. 


12T2  -t-  13T3) 


dT- 

dz 


=  B'  ( 


T 


2  j  -  5 
+  -  T. 


6 j  -  3  ^"2  2j  -  1  ^_1 


4  (  j  -  1 )  6j+l 

-  T  +  - 


2j  -  1 


6j  -  3 


Tj+1> 


dT 


10 


dz 


=  6  ((97  -  102)  Tg  -  (507  -  1060)  T( 


where 


+  (2257  -  2910)  T1q  +  y&Tg) 


B  *  =  B/Ar‘ 


P 


60  Ar  hf 


k, 


7 


1888/(184  +  p) 
B'/570 


r=R 


(III. 7) 


6 
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to  attempt  by  hand  calculations,  is  now  possible  and  practic¬ 
able  with  digital  or  analog  computers. 

The  problem  is  treated  as  an  initial  value  problem, 
which  would  be  the  case  for  cocurrent  flow.  All  the  variables 
are  known  at  some  one  value  of  the  independent  variable,  namely 
at  the  reactor  inlet,  i.e.  z  =  0.  For  countercurrent  reactors 
the  initial  values  of  X,  P  and  Tg  would  have  to  be  assumed  and 
the  solution  is  obtained  by  trial-and-error .  The  Runge-Kutta- 
Gill  method,  outlined  in  Appendix  B,  is  used  for  the  numerical 
integration . 

Before  numerical  examples  are  presented,  an  alter¬ 
native  approach  will  be  discussed.  One  of  the  most  powerful 
tools  for  solving  sets  of  simultaneous  ordinary  differential 
equations  is  the  analog  computer,  which  by  its  inherently 
parallel  operation  avoids  the  problem  of  stability  encountered 
in  finite-difference  methods.  The  linear  differential  equations 
can  be  simulated  on  an  analog  computer  with  a  high  degree  of 
accuracy  and  require  comparatively  little  equipment  -  about 
four  amplifiers  per  equation.  The  non-linear  equations  require 
the  use  of  special  function  generators  which  have  a  lower 
degree  of  accuracy  and  require  a  great  many  more  electronic 
components.  In  fact,  it  is  fair  to  say  that  the  difficulty 
of  generating  non-linear  terms  is  the  biggest  physical  problem 
in  using  analog  computers.  There  is,  however,  no  inherent 
theoretical  limitation  to  the  solution  of  such  sets  of 
differential  equations  by  analog  computers. 
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As  a  check  on  the  accuracy  of  the  methods  proposed, 
two  numerical  example  problems  have  been  solved. 

III.l  Example  Problem  No,  1 

This  example  problem  on  moving  bed  heat  exchanger 
was  presented  in  Amundson's  paper (10): 

Limestone  is  calcined  in  a  continuous  countercurrent 
vertical  lime  kiln.  The  calcium  carbonate  is  fed  in  at  the 
top  in  2"  pieces  at  a  superficial  mass  velocity  of  230  pounds 
per  square  foot  per  hour.  Flue  gas  and  carbon  dioxide  from 
the  calcination  rise  through  the  bed  at  a  superficial  mass 
velocity  of  2500  pounds  per  square  foot  per  hour.  If  the 
limestone  enters  at  100 °F  and  the  gas  leaves  at  200 °F,  what 
is  the  temperature  distribution  in  the  kiln  as  a  function  of 
distance  from  the  top,  assuming  no  chemical  reaction? 


Physical  data  are 
specific  heat  of  gas 
specific  heat  of  calcium  carbonate 
density  of  calcium  carbonate 
heat  transfer  coefficient 

thermal  conductivity  of  calcium  carbonate 
fractional  void  space 

Since  there  is  no  chemical  reaction  present  and 
the  kiln  is  assumed  to  be  adiabatic,  without  consideration  of 
the  pressure  drop,  the  gas  phase  is,  therefore,  described  as 


=  0.25  BTU/ ( lb) ( °F ) 

=  0.28  BTU/ ( lb) ( °F ) 

=  162  lb/ (ft3) 

=  78  BTU/(hr) (ft2) ( °F) 

=  1.3  BTU/(hr)  ( ft  )  (  °F) 

=  0.50 
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3(1  - e  )  hf  (Tg  -  Ts) 
G  cf  R 


(III. 8) 


Together  with  equation  (III. 7),  they  form  a  set  of  simultaneous 
linear  first  order  ordinary  differential  equations.  A  Pace 
231R  analog  computer  was  used  to  obtain  solution.  The  circuit 
diagram  for  this  case  is  shown  in  Figure  2  and  the  result 
generated  by  the  X-Y  plotter  is  given  in  Graph  No.  2.  The 
numerical  solution,  obtained  by  the  Runge-Kutta-Gill  method, 
was  performed  on  an  I.B.M.  1620  digital  computer.  The  pro¬ 
gram,  written  in  Fortran  IV,  is  listed  in  Appendix  F. 

This  problem  has  been  solved  by  Furnas (11)  analyti¬ 
cally  in  neglecting  conduction  in  the  solid,  by  Lovell  and 
Karnof sky ( 12 )  using  a  modified  Schmidt  method,  and  by 
Amundson(lO)  using  successive  Laplace  transformations  to 
give  the  exact  solution.  Results  from  the  above  methods  are 
plotted  in  Graph  No.  3  for  comparison.  It  is  obvious  that  the 
results  of  the  proposed  method  very  nearly  coincide  with 
those  by  Amundson.  As  a  better  comparison  the  results  are 
tabulated  in  Table  2.  This  confirms  the  accuracy,  both  of 
the  finite  difference  expressions  and  of  the  integrative 
technique  used. 


I  ilv  bcus  £  dnw^f  ‘2  ni  nwoxfa  ai.  qbbo  ai/ii  “tc  -  • r- 1  '•* 
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Figure  2 

Analog Computer  Circuit  Diagram  f or 
Moving  Bed  Heat  Exchanger. 
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Graph  No.  2 

Temjp^ratu re  Profile  for  Moving  Bed  Heat 
Exchanger  from  Analog  Computer. 


Equivalent  Length  of  Heat  Exchanger  xlOO 
(K.t/R2)  or  (K.z/R2.Us) 


Reduced  Temperature  of  Gas  (Tg-T0)/(Tg- 


Graph  No.  3 
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Temperature  Profiles  along  the  Moving-bed  Reactor 

in  countercurrent  system 


with  (Gs. Cs/Ga. Cf )  =  1.0 
(ks/hf.R)  =0.2 


Equivalent  Length  of  Reactor 

(K.t/R2)  or  (K.z/R2.Us) 
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TABLE  2 

Comparison  of  Solutions  in  Moving- Bed 

Heat  Exchanger  Example 


Equivalent 

Length 

of 

Exchanger 

(Kz/R2uJ 


T  -  T 
ig  -Lo 

Reduced  Temperature  (— - — ) 

1 go  ~  1 o 

Analytical  Explicit  Explicit 

Solution  by  Solution  by  Solution  by 

Amundson  Digital  Computer  Analog  Computer 


0.05 

1.46 

1.458 

1.47 

0.10 

1.85 

1.848 

1.85 

0.15 

2  .23 

2  .227 

2.22 

0.1854 

2.50 

2.494 

2.500 

IBM  7040 

Pace  231R 

Fortran  IV 

20  sec.  for 
63  increments 
in  z 
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III .2 _ Example  Problem  No .  2 

This  problem,  taken  from  Walas '  "Reaction  Kinetics 
for  Chemical  Engineers" ( 13)  is  given  as  follows: 

Butane  is  cracked  in  a  moving  bed  pebble  heater, 
according  to  the  following  reaction, 

C4H10  - *  2C2H4  +  H2  (III. 9) 

with 

AH^Qop  =  98,600  BTU/lb-mole  converted 
AHi2oq°f  ”  99,100  BTU/lb-mole  converted 
The  reaction  is  first-order  with 

log  k  =  -22,100/T  +  12.45 

v/here  is  in  °R  and  k  in  sec"1'. 

The  pebbles  are  alumina  spheres  5/16  inch  in  diameter  and 
have  a  heat  capacity  of  25  BTU/(cu . ft . ) ( °F) ,  a  surface  area 
of  126  sq.ft ./(cu.ft .) ,  and  a  void  fraction  of  0.45.  After 
preheated  uniformly  to  2000 °F  the  pebbles  are  allowed  to  flow 
downward  through  the  reactor  in  parallel  to  butane,  which 
enters  at  500°F  and  20  psig.  The  butane  flow  rate  is  42.6 
lb-moles/(hr ) ( sq. ft . )  and  the  pebble  rate  is  682  cu.ft.  per 
(hr) (sq.ft .) .  The  heat  transfer  coefficient  from  the  gas  to 
the  pebbles  is  97.6  BTU/(hr ) ( sq . ft . ) ( °F) .  It  is  desired  to 
calculate  the  length  of  reactor  required  to  obtain  various 
fraction  conversion  of  the  butane.  Property  values  of  the  gas 


are  given  in  the  text . 
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For  the  first-order  reaction  (III. 9)  if  the  gas  is 
assumed  to  behave  ideally,  then  the  rate  of  reaction  is  given 
by 


where 


R  (x,p,Tg)  =  kCa 


k  P  (1  -  X) 
R'  Tg  (1  +  2X) 


lb-mole 
^hr  ft^  ^ 


(III. 10) 


Ca  is  the  concentration  of  n-butane 
X  is  the  fractional  conversion  of  n-butane  in 
terms  of  total  feed 
P  is  the  total  pressure 


dx  6  P  •  ( 1  -  X) • exp  (12.45  -  22100/T_) 

—  _  s? 

dz  Gm°'R'  Tg ' ( 1  +  2X) 


(III. 11) 


(HI.l) 


If  Tb  is  the  reference  temperature  then  equation 
can  be  written  as 

hf  As  f 

- —  (Ts  -  Tg)  -  (-CPj^  +  2Cp2  +  Cp3)  (Tg-Tb)+AHTb- 

Gm  * 


dx 

dz 


dz 


((1  -  X)  Cp1  +  2XCp2  +  XCp3) 


(III. 12) 


where  Cp^,  Cp2,  C^3  are  ^ie  mo-ar  specific  heats  of  n-butane, 
ethylene  and  hydrogen  respectively. 

The  density  at  any  point  in  the  reactor  for  ideal  gas 


behavior  is  given  as 
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58 

Pf  T a  14~  7 

(1  +  2X)  (379)  (-— )  ( - ) 

520  P 


5.4P 


(1  H*  2X)  Tg 


(III. 13) 


so  that  with  mean  values  for  \±  equation  (III.  3)  can  be 
simplified  into 

dP  (1  +  2X)  Tq 

—  -  0.01  x  - - —  —  (III.  14) 

dz  p 


As  the  solid  radius  is  not  too  large  it  is  sufficient 
to  divide  it  into  5  equal  sections  so  that  equation  (III.  5) 
is  expressed  into  a  set  of  but  5  equations,  listed  in  Table  3. 

The  complete  set  of  equations  for  computation 
includes  (III. 15)  for  the  solid  and  (III. 11),  (III. 12), 

(III. 14)  for  the  gas.  It  can  be  seen  that  the  equations 
describing  the  gas  phase  are  non-linear,  so  that  analytical 
solution  is  impossible.  However,  in  the  solution  by  Walas, 
the  internal  temperature  of  the  pebbles  is  considered  to  be 
uniform  and,  therefore,  the  thermal  conductivity  of  the 
pebbles  does  not  enter  into  the  calculations.  Thus,  with 
finite  increments  along  the  reactor  length,  a  solution  by 
trial-and-error  for  the  gas  temperature  at  the  end  of  each 
increment  can  be  obtained.  In  the  present  proposed  metlioc  , 
the  set  of  equations  for  the  system  wa 


s  solved  numerically  by 
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r=0 


where 


TABLE  3 


The  Computational  System  for  Solid  in  the 

Moving  Bed  Reactor 


r=R 


dT-,  B' 

— -  = - (-26TX  +  27T2  -  T3) 

dz  8 


dT9  B' 

— -  =  —  ( -T-i  -  12T7  +  13To) 

dz  9  1  ^ 


dTo  B' 

-  =  -  (2Tn  +  3T0  -  24T»  +  19TJ 

dz  15  1  1  5  * 


dT,  B* 

~  =  —  (2T9  +  9T3  -  36T4  +  25T5) 

dz  21 


dTc 

— —  =  6  ((9y  -  42)T3  -  (507  -  460)T4 

dz 

+  (2257  -  1410)T5  +  7^Tg) 

B'  =  B/Ar2 

p  =  (60  Ar  hf)As 

7  =  983/(184  +  p) 

6  =  B 1 /( 270) 


(III. 15) 
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the  Runge-Kutta-Gill  method,  using  a  pebble  thermal  con¬ 
ductivity  of  1.8  BTU/(hr . ft . °F) . 

The  results  of  the  different  methods  are  illustrated 
in  Graphs  No.  4a  and  4b  for  the  fractional  conversion  and  the 
gas  temperature  profiles  respectively.  The  differences  are 
considerable.  By  comparing  the  results  tabulated  in  Table  4a, 
from  Walas,  with  those  in  Table  4b  from  proposed  explicit  method, 
the  difference  becomes  obvious.  The  temperature  difference 
between  gas  and  solid  is  as  high  as  300 °F  for  the  greater 
part  of  the  reactor  in  the  former,  while  in  the  latter,  it 
is  less  than  100 °F.  In  fact  the  temperature  gradients  in  the 
solid  can  be  found  to  be  at  least  15 °F  for  half  of  the  reactor 
length,  which  is  certainly  not  negligible. 

The  program  was  written  in  Fortran  IV  and  listed 
in  Appendix  F. 


III. 3  Summary 

Accurate  and  stable  methods  have  been  proposed 
for  the  solution  of  one-dimensional  moving  bed  reactors, 
where  chemical  reaction,  if  any,  occurs  only  in  the  gas 
phase.  Finite-difference  techniques  lead  to  a  set  of 
simultaneous  first-order  ordinary  differential  equations, 
which  may  then  be  solved  either  by  using  the  Runge-Kutta- 
Gill  method  on  a  digital  computer  or  by  an  analog  computer. 


Fractional  Conversion 
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Graph  4a 

Fractional  Conver s ion  Profiles  in 
Thermal.  Cracking  of  n- Butane  in  Moving  Bed _ Reactor. 


Length  of  Reactor  z  (ft) . 


Gas  Temperature 
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Graph  4b 


Gas  Temperature  Profile  in 

Thermal  Cracking  of  n-Butane. 


Length  of  Reactor  z  (ft)  . 
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TABLE  4a 


Results  in  Thermal  Cracking  of 

Butane  from  W a 1 a s 


Reactor 

Length 

(ft) 

Gas 

Temp . 

(°R) 

Solid 
Temp . 

(eR) 

Temp . 

Difference 

(  °R) 

Conversion 

0 

500 

2000 

— 

0 

0.105 

1200 

1919 

- 

0 

0.17  5 

1340 

1901 

561 

0.003 

0.212 

1408 

1889 

481 

0.015 

0.259 

1464 

1874 

410 

0.05 

0.305 

147  3 

1864 

391 

o 

* 

t— 1 

o 

0.520 

1496 

1809 

313 

0.30 

0.792 

1520 

17  56 

234 

0.50 

1.179 

1541 

1703 

162 

0.70 

1.988 

1575 

1650 

75 

0.90 
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TABLE  4b 

Results  in  Thermal  Cracking  of 

Butane  by  Proposed  Method 


Reactor 

Length 

(ft) 

Gas 

Temp . 

(  °R) 

Surface 
Temp . 

(  ®R) 

Temp  . 

Difference 
( °R) 

Conversion 

0 

500 

2000 

- 

0 

0.10 

1084.4 

1768.9 

- 

0 

0.30 

1432.7 

1713.1 

280.4 

0 • 0663 

0.60 

1451.4 

1640,2 

188.8 

0.2582 

1.20 

1452  .9 

1547  .8 

94.9 

0.4941 

1.80 

1444.7 

1497  .7 

53.0 

0.6229 

2.40 

1435.0 

1467  .7 

32.7 

0.7008 

3.00 

1426.2 

1448.1 

21.9 

0.7521 

3.60 

1418.7 

1434.4 

15.7 

0.7883 

4.20 

1412.4 

1424.2 

11.8 

0.8152 

4.80 

1407.2 

1416.4 

9.2 

0.8359 
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Both  example  problems  demonstrate  the  validity  and 
practicability  of  the  proposed  numerical  method.  From  the 
object  times  of  computation  the  method  can  be  considered  as 
efficient.  The  second  problem,  however,  reveals  the  error 
of  assuming  uniform  temperature  gradients  in  the  solid, 
which  leads  to  considerable  underdesign  of  the  reactor. 
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IV.  FIXED  BED  REGENERATIVE  REACTOR  DESIGN 

Under  the  assumptions  outlined  in  the  introduction, 
the  reactor  may  be  described  mathematically  as  follows: 

1.  The  energy  balance  on  the  fluid  contained  in  a 
differential  volume  of  the  reactor  of  unit  cross-sectional  area 
and  length  dz  is  given  by 


d 

£<32  —  (Tg  Cf  Pf  vf)  +  As  dz  hf  (Tg  -  Ts) 


N 


+  G, 


m 


o 


•  dz  • 


l 


(IV.l) 


which  simplifies  to 


G  —  ( c  f  T  ) 
oz 


3(1  -  €)  hf  (Tg.  -  Ts) 


R 


N 


+  G, 


m 


o 


I 


(IV. 2) 


where 


x j  ( X-j^ ,  X2  / 


•  •  • 


j—1,2, ••• 


N 


(IV. 3) 
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Equation  (IV. 3)  stands  for  the  fractional  conversion  rate  of 
the  jth  limiting  component  in  the  system  of  N  independent 
reactions . 

a 

2.  The  pressure  drop  for  fixed  bed,  if  considered, 
is  given  empirically ( 14)  as 


dp 

dz 


-G(  1  -  £)  ^  150  (1  -  Op. 

pf  Dp  €3  Dp 


+  1.7  5G) 


(IV.  4) 


3.  The  thermal  behavior  of  the  individual  solid 


particle  is  described  by 


c  s  P  s 


dT 

dt 


~v2  „ 

o  T 


ho  ( 


dr4 


+ 


2  dT 
r  dr 


(IV.  5) 


Again  equations  (IV. 2)  and  (IV. 5)  are  coupled  by  a  heat 
balance  at  the  gas-solid  interface. 


dT 


ks  ^  dr^  r=R 


hf  (Ts  -  T  ) 


The  other  boundary  condition  for  the  solid  phase  is 


OT 


dr 


=  0 


at  r  =  0 


which  denotes  radial  symmetry. 

The  initial  conditions  are 


Tg  =  Tg 


T  =  T  (r ,  z) 


at  z  =  0 


at  t 


dz 

v^ 


0 


(IV. 6) 


(IV.  7) 


(IV. 8) 
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Equation  (IV, 2)  is  based  on  a  heat  balance  on  a 
fixed  element  in  space,  i.e.  Eulerian  view  point.  Alternately, 
by  adopting  the  more  convenient  Lagrangian  point  of  view,  that 
is  to  follow  a  fixed  differential  mass  of  fluid  as  it  proceeds 
through  the  reactor,  it  is  possible  to  eliminate  the  time 
variable  and  reduce  equation  (IV. 2)  to  an  ordinary  differential 
equation.  This  implies  that  there  is  no  change  in  mass  of 
the  fluid  in  passing  through  the  differential  element  of 
reactor  with  time.  In  terms  of  the  new  transformed  variable 
of 


the  equations  describing  the  system  are  now 


dT, 


dz 


3(1  -  «  )hf 


R  G  c 


„  °  N  , 

Gm  v-1  dxi 

s  *  Gc.  L  dz  3 


(IV.  10) 


j=l 


together  with  equations  (IV. 3)  and  (IV. 4)  for  the  gas  phase; 


dT 

dQ 


=  K  (- 


d2T 


2  dt 

+  —  — ) 
r  dr 


(IV. 11) 


and  the  boundary  conditions  (IV. 6)  and  (IV. 7)  for  the  solid. 
While  the  initial  conditions  are 
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and 


rp 

"g 


(t) 


j  j 

p  =  p° 


T  =  T  (r ,  0) 


at  z  =  0 


9  =  0 


(IV. 12) 


IV .I  Previous  Analytical  Solutions 


A  general  analytical  solution  for  the  above 
coupled  system  of  differential  equations  is  not  available. 
However,  an  analytical  solution  of  a  simplified  case  was 
proposed  by  Schumann (15) .  In  the  absence  of  chemical  reaction 
in  the  gas  phase  and  neglecting  the  temperature  gradient  in 
the  solid  the  mathematical  model  can  be  reduced  to  the  fol¬ 
lowing 

oTq  cTq  -3(1  -  € )  hf 

—  +  v f  —  =  — - - (T  -  T)  (IV. 13) 

St  Sz  G,  R  cf  Pf 


Sr 

St 


R  c 

s 


T) 


(IV.  14) 


In  terms  of  the  new  variables,  defined  below, 


Y  = 


3  (1  -  €  )  h£  z 

G  R  c£  pf  v£ 


and 


Z 


3h£ 

^  cs  Ps 


(t 


(IV. 15) 
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equations  (IV.  13)  and  (IV. 14)  become 


dT, 


<3y 


=  T  -  T, 


(IV. 16) 


T  -  T 

g 


(IV.  17) 


with  initial  conditions 


T  -  T 

g  g 


T  =  T 


o 


at  Y  =  0 

at  Z  =  0 


(IV.  18) 


Strictly  speaking,  the  above  equations  are  true 
only  if  Vf  is  constant,  but  Jakob(16)  has  shown  that  for  most 
practical  cases,  the  error  introduced  by  moderate  variations 
in  V  -  is  negligible. 

After  solving  equations  (IV.  16)  and  (IV.  17)  simul¬ 
taneously  with  the  initial  conditions  (IV. 18),  the  solution 


for  the  gas  temperature  is 

rp  _  mO 

_ 


OO 


-(Y+Z) 


n 


T  °  -  T° 

g 


Z  Mn  (YZ) 


(IV.  19) 


n=0 


and  for  the  solid  temperature  is 


T  -  T 


T  °  -  T° 

g 


=  e 


-(Y+Z) 


oo 


n=l 


Zn  M  (YZ) 
n 


(IV. 20) 


In  letting  a  =  ZY,  then 
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2 


3 


a 


a 


Mo(a) 


1  +  a  + 


+ 


+ 


oo 


(IV. 21) 


and 


da 


IV. 2  Physical  Representation  of  Mathematical  Model 

The  mathematical  model  proposed  for  the  numerical 
analysis  consists  of 

1.  an  elemental  volume  of  reactor  of  unit  cross- 
section  and  length  Az,  in  which  the  temperature  distribution 
in  the  solid  particle  is  the  same  throughout  the  element. 

2 .  an  elemental  slug  of  fluid  of  such  size  that 
it  would  pass  through  the  reactor  element  in  time  A9 . 

The  model  merely  involves  bringing  the  solid  and 
the  fluid  into  contact  at  time  zero  and  following  the  behavior 
of  both  over  the  time  interval  AG. 

The  slug  of  fluid  then  goes  on  to  contact  the  next 
reactor  element  and  so  on  through  the  reactor.  In  this  way 
the  thermal  condition  of  the  solid  bed  is  determined  for  one 
AG.  The  procedure  is  then  repeated  for  the  next  slug  of  fluid. 
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Figure  _3 


Diagrammati 


the  MathemaL ti c a  1  Model  for  the 


Regenerative  Bed  Reactor, 


Gas  element 


Intersections  may  be  considered  to  be 
stirred  tanks. 
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It  should  be  noted  that  in  the  following  calculations,  the 
time  average  is  used  for  the  solid  surface  temperature  and 
position  average  for  the  gas  temperature.  In  fact,  the  model 
reduces  to  a  series  of  stirred  tank  calculations.  Diagram- 
metically  it  is  illustrated  in  Figure  3. 

IV. 3  The  Fixed  Bed  Regenerative  Heat  Exchanger 

In  this  case,  there  is  no  chemical  reaction  and  if 
the  pressure  drop  is  not  to  be  considered  then  equation  (IV. 10) 
together  with  equation  (IV. 11)  form  a  linear  system.  That 
means  a  completely  stable,  implicit  method,  such  as  outlined 
in  Section  II,  can  be  used  for  the  numerical  integration. 


The  finite-difference  approximation  of  equation 


(IV. 10)  is  given  by 


Az  •  3(1  -  )hf 


( T  -  T 

sj+h,i+h  9±+h 


(IV. 22) 


\ 


where 


Az  3(1  -  £)hf 


(IV. 23) 


A 


1 


R  G  c 


f 


In  the  above  equations  the  subscript  j  is  used  to 
denote  reactor  position  and  the  subscript  i  to  denote  time. 


1-41 


Thus,  the  solid  surface  temperature  is  at  the  time  average 
and  the  gas  temperature  at  the  position  average  as  described 
in  the  model. 

For  the  solid  phase,  the  partial  differential 
equation  (IV. 11)  is  replaced  by  a  set  of  ordinary  differential 
equations,  listed  in  Table  3,  with  the  radius  being  divided 
into  5  equal  lengths  and  the  temperature  at  the  center  of 
each  length  representing  the  mean  temperature  of  the 
spherical  shell.  Expressing  equations  (IV. 22)  and  (IV. 11) 
in  Crank-Nicholoson  implicit  form,  the  corrputational  system 
becomes: 


13B  27B  B 

^  +  -r!Ti,i+i  -  <ir)T2,i+i  +  (ii)T3,i+i 


13B  27B  B 

=  (1  -  - )T,  .  +  (— -)T0  .  -  (— )T  . 

8  1/i  16  2 ,  i  16  3,i 


4)Ti.i+i  +  (1  +  r)T2-i+i 


2B 


13B 

(~)T3,i+l 


B 


2B 


13B 


=  -  (— )Tn  .  +  (1  -  - )T  .  +  (— )T  • 

18  3  ^ ,  l  8  3,i 


B  B  4B 

( - )  T  ,  -  ( - )T0  .  _  +  (It  — )T 

15  1,1+1  10  2,i+l 


5  3,i+l 


19B 

( - )  T  . 

30  4,i+l 


B  B  4B  19B 

(T7)Tl,i  +  (V)T2,i  +  (1T)T3-i  + 
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B 

*2l*T2,i+l 


3B 

<14>T3,i+l 


GB 

(i  +  y :i 


.>  sr> 

(  )  T 

'  <!  2  ;  5,  i-M 


,B 

( — )rJ7  . 

71  2,:i. 


/3Bx 

(  — )T- 

]  A.  3  f  1 


(1  '  7-)T4,x 


■  2  5B 

+  (~ . )Tr  . 

42  -1- 


540 


(42  -  9y)  T-j  -I-  (50y  -  460 )T.  .  fl  f  ( - -V  1410  -  2257)Tr  ._L. 

1  ft,  1*1* X  g  Jt  .1  +  1 


-  7pVi 


(IV. 24) 


540 

(97  -  42)T  .  +  (460  -  507)17  .  -I-  ( -  -  1410  *1-  2257)17  . 

*  •  1  1  b  1 


+  7PT 


-  9T 


3,i+l 


SOT 


4,  i+1 


225T 


5/i+l 


(184  +  2(184  -I-  p)/A3)T 


q  • 


]+l 


=  9T 


3#  i 


50T 


4,i 


225T 


5,i 


(184  -  2(184  +  p)/A,  )Tct 

1  9j 


where 

Ar  =  R/5 

B  =  A©K/Ar 2 

p  =  60ArhfAs 


7 


992/(184  +  p) 
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It  should  be  noted  that  all  the  terms  in  equation 
(III. 11)  are  third  order  correct.  A  system  of  second  order 
correct  expressions  had  been  tried  but  this  led  to  unsatisfac¬ 
tory  results.  The  system  of  equations  (IV. 24)  is,  in  fact,  a 
set  of  simultaneous  linear  algebraic  equations  which  can  be 
written  algebraically  as  in  Table  5.  The  coefficient  matrix 

t 

can  be  transformed  by  simple  arithmetic  operations  into  a 
tridiagonal  matrix,  which  can  then  be  readily  solved  by  the 
Thomas  method (17) . 

A  numerical  problem  was  solved  to  check  the  feas¬ 
ibility  of  the  implicit  method  against  the  analytical  solution 
by  Schumann  for  a  system  with  negligible  intraparticle 
temperature  gradients. 

IV o 4  The  Regenerative  Bed  Fluid  Phase  Reactor 

The  presence  of  non-linear  terms,  as  a  result  of 
chemical  reaction  in  equation  (IV. 10),  make  difficult  if  not 
outright  impractical  the  implicit  method  outlined  in  the  pre¬ 
ceding  section.  Thus,  resource  to  some  kind  of  explicit 
method  operating  on  equation  (IV. 10)  seems  loqical.  However, 
the  stability  of  these  methods  depends  upon  the  magnitude 
of  the  term  (KA9/Ar2),  which  cannot  exceed  some  limiting  value. 
This  value  cannot  be  easily  predicted  for  systems  of  non-linear 
differential  equations  and  depends  on  the  very  nature  of  the 
non-linear  terms.  Usually,  in  order  to  preserve  stability,  AO 
has  to  be  made  so  small  that  an  extraordinarily  large  number 
of  time  steps  must  be  used  to  obtain  a  practical  solution. 


TABLE  5 


The  Simultaneous  Linear  Algebraic  Form  of  Equations 

for  Fixed  Bed  Regenerative  Heat  Exchanger 


—  — 

all  a12  a13 

Tl, i+1 

a21  a22  a23 

T2, i+1 

a31  a32  a33  a34 

T3, i+1 

a42  a43  a44  a45 

T4, i+1 

a53  a54  a55  a56 

T5,i+1 

a63  a64  a65  a66 

T6, i+1 

_  _ 

—  — 

where 


T 

■fa, i+1 


T 


2,  i' 


k 


1,  2, 


•  *  • 


6 


■ 
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It  is  proposed  here  that  the  problem  can  be  reduced 
to  manageable  proportions  by  decoupling  the  two  phases.  In 
principle,  the  method  involves  substituting  an  explicit  function 
for  Ts  in  terms  of  9  in  place  of  the  coupling  boundary  condition 
given  by  equation  (IV. 6).  As  a  first  approximation,  it  is 
assumed  that  over  a  short  time  increment  the  surface  tempera¬ 
ture  changes  linearly  with  time,  i.e. 

Ts  =  Ts°  -I-  q  Q  (IV. 2 5) 

where  the  superscript  "o"  signifies  the  beginning  of  the  time 
increment  and  q  is  some  constant  to  be  determined.  According 
to  the  model,  the  elemental  slug  of  fluid  will  see  a  constant 
solid  surface  temperature  of  (T  °  +  qAQ/2). 

O 

Thus,  the  correct  value  of  q  is  that  which  will 
equate  Qs  the  total  heat  gained  by  the  solid  phase  in  time  A© 
and  Qg  the  enthalpy  lost  by  the  elemental  slug  of  fluid  under 
consideration.  To  evaluate  these  quantities,  equations 
( IV  e  3 ) ,  (IV. 4),  (IV.  10)  and  (IV. 11)  have  to  be  solved  for 
the  assumed  boundary  condition . 

The  solid  phase  will  be  considered  first.  Since 
equation  (IV. 11)  is  linear,  the  principle  of  superposition 
can  be  applied;  hence,  it  can  be  replaced  by  two  other 
systems,  each  with  simpler  boundary  conditions.  Let. 

T=u+w  (IV. 26) 


=  K  ( 


d2w 

aP 


2  dw 
+  -  — ) 

r  or 


Then 


dw 

d© 


(IV. 27) 


with  boundary  conditions 


dw 

0 

—  = 

at  r  =  0 

( XV  ,28a) 

dr 

o 

(IV, 28b) 

w  = 

Ts 

+ 

q  © 

at  r  =  R 

and  initial 

condition 

w  = 

0 

9 

=  0 

(IV. 28c) 

Also 

du 

>2 
o  u 

2 

du 

—  = 

K 

(— T 

+  — 

— ) 

(IV. 29) 

do 

Sr2 

r 

dr 

with  boundary  conditions 


W\¥ 

II 

0 

n 

ii 

o 

(IV.  30a) 

U  = 

0 

r  =  R 

(IV.  30b) 

and  initial  condition 

u  = 

T(r) 

o 

ii 

a> 

(IV. 30c) 

Equation  ( IV 

.27)  with  its 

associated  boundary  con- 

ditions  can  readily  be  solved  by  applying  the  Laplace  trans¬ 
formation  to  obtain  a  series  solution  involving  the  error 
function  as  follows: 


RT, 


OQ 


n, 


erfc  ( — ) 
* 


W 


r 
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+ 


oo 


n=0 


2K 


Y 

erfc  (— ) 
A 


( — )  ^  exp 
7TK 


II 

A 


(9  + 


2K 


erfc 


9  3s 

6  ( — )  exp 

7TK 


(-( 


A 


(IV. 31) 


where 

Yn  =  (2n  +1)  R  -  r 

A  =  2  (K9)^,/2 

8r  =  (2n  +1)  R  +  r 

Equation  (lV.29)  with  its  associated  boundary  condi¬ 
tions  is  usually  handled  by  the  method  of  separation  of 
variables  with  the  initial  condition  in  discrete  form. 
Expansion  into  an  infinite  series  of  orthogonal  functions, 
results  in 


R 


u  = 


2 _ 

Rr 


(IV. 32) 

y  exp  (-Kan2©)-  sin  anr  •  J  r'T(r')-  sin  anr'-  dr' 


n=l 


where 

a  =  mr/R 

n 

However,  practical  computational  problems  arise 
since  the  series  which  gives  the  heat  content  change  over  a 
small  interval  of  time  A9,  as  shown  below. 
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A© 


du 

U~)  0 

or  r=R 


d© 


oo 


(-l)n 


V 


2  2  A  _ 
Kn  7T  A© 

exp  ( - - ) 


R' 


r'T(r') 


nTTr  * 

sin  - —  dr  J 

R 


(IV. 33) 


converges  very  slowly  due  to  the  smallness  of  the  quantity 
(Kn27T2A0/R2)  . 


To  resolve  this  difficulty,  Green's  function  for  an 
instantaneous  spherical  surface  source  is  used.  For  an  infinite 
sphere  the  temperature  at  r,  due  to  an  instantaneous  spherical 
surface  source  of  strength  Q'  at  ©  =  0  and  of  radius  r’,  is 
given(18)  as 


u 


1 


Q’ 


8Trrr  5  (ttKQ) 


172 


- (r  -  r  ’ ) 2 

<  exp  ( - — ) 

4K© 


-(r  +  r • ) 2 

exp  (— - - - ) 

4K© 

(IV. 34) 


Now  for  a  finite  sphere  with  unit  instantaneous 
spherical  surface  source  at  ©  =  0  and  of  radius  r •  the 
temperature  is  given  by  u 1 ,  such  that 


u*  =  ux  +  u2  (IV. 35) 

where  u2  vanished  at  ©  =  0  and  is  such  that  u'  satisfies 

du '  d2  u '  2  du ' 

-  =  K  (— —  +  ~  - ) 

air2 


3© 


r  Sr 


(IV.  36) 
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with 


u  ’ 

u ' 

du 

dr 


=  unit  inst.  source 


=  0 


=  0 


9  =  0 
r  =  R 

r  =  0 


Since  U2  also  satisfies  the  heat  conduction  equation 


d  (ru2) 
dr2 


1  d(ru2) 
K  d9 


0 


(IV.  37) 


Taking  the  Laplace  transformation  with  respect  to  9,  (IV. 37) 
becomes 


d2(ru2)  _ 

-  -  qz  (ruJ 

dr2 


0 


( TV „ 38) 


where 


oo 


u. 


e”^9  u2  d9 


0 


q2  =  p/k  . 


The  general  solution  of  equation  (IV. 38)  is 


u- 


(c-^  sinh  qr  +  c2  cosh  qr)  #  (IV. 39) 


At  r  =  0,  u2  is  finite,  so  c2  =  0 


However,  at  r  =  R 


u*  =  ui  +  u2 


0 


.  > 
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Therefore , 


u 


1 


1 


8n"Rr '  Kq 


exp 


(-q(R 


r*)) 


exp  ( -q(R  +  r ' ) )  » 


+ 


c^(sinh  qR)/R 


=  0 


where  Q'  is  unity. 

After  solving  for  c-^  and  back  substituting  into  u ' , 

it  gives 


sinh  qr*sinh  q(R  -  r') 

u '  =  -  (IV. 40) 

47rrr'Kq  sinh  qR 


Expanding  in  a  series  of  negative  exponentials  and 
taking  the  back  transformation  from  tables,  the  solution  is 


u '  = 


8rr  'ttCttK©) 


1/2 


oo 


n=-oo 


exp 


( 2nR  -  r  +  r')2 


4K9 


(IV. 41) 

> 

If  now  the  initial  condition  for  u  is  considered, 

2 

i.e.  of  source  strength  4TT(r')  T(r')dr’  on  the  sphere  at  r' 


exp 


( 2nR  +  r  +  r')2 


4K9 
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oo  R 
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1  / 

- — ^  )  f  r-T(r'). 
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n=-oo  0 
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4K0 


exp 


( 2nR  +  r  +  r 1 )  ^ 


4K0 


*  dr 


(IV. 42) 


In  essence,  the  Green's  function  permits  the 
evaluation  of  the  temperature  at  any  point  due  to  an  instant¬ 
aneous  surface  source  which  moves  from  the  centre  to  the  outer 
boundary o  Mathematically  this  is  equivalent  to  calculating 
the  additive  component  of  temperature  due  to  the  initial 
temperature  distribution  within  the  particle. 

Thus,  the  solid  temperature  T  is  obtained  simply 
by  adding  equations  (IV. 31)  and  (IV. 42),  that  is 
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,  Yn  Yn  9  ^  Yn  2 

(0  +  — = — )  erfc  (— )  -  y  ( — -)  2  exp  (-( — -)  ) 
2K  A  n  7TK  A 

(IV.  43) 

where 

Yn  =  (2n  +  1)  R  -r  6n  =  (2n  +  1)R  +  r 

A  =  2(K0)^ 

It  is  not  difficult,  though  tedious,  to  verify  that  this 
solution  indeed  satisfies  the  heat  conduction  equation  and 
the  associated  boundary  condition0.  However,  it  is  proved  in 
Appendix  E,  that  it  satisfies  the  initial  condition  T(r) . 

Now  the  total  heat  content  change  Qs  of  the  solid 
over  the  time  interval  A0  can  be  calculated  by  integrating 
the  temperature  distribution  thus  obtained  over  the  whole 
volume  of  the  particle.  However,  this  would  have  to  be 
repeated  for  each  new  value  of  q. 

A  more  efficient  method  is  to  evaluate  the  heat  flux 
at  the  surface  and  integrate  this  over  the  time  interval. 


Q 


s 


3(1  -  t )  Az 


R 


dT 


d9 
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3(1  -  OAz  ks 
— - —  < 

R 


oo 


n=0 


r'T(r') 
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2  /kA© 


2nR  +  R  -  r' 

-  erfc  ( - -  -  - - ) 

2  J  KA  © 
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+  q 


-A©2 

2R 


+ 


4 

3 


+ 


n2R2 

exp  ( - ) 

KA© 


2nR 

erfc 

K 


nR 


'  /ka©} 


n=l 


KA©n 

— 7T-  erfc 
2R2 


nR 


(-7=) 

7ka© 


-n2R2 


A©  K  2 

—  —  n  exp  ( - 

Kit  R  KA© 


)  +  erfc  ( — --  — ) 

Vka© 


nR 


16q  R 


3  00 


3K  ^  n=l 


n 


nR 


erfc  (- - )  •  pf 


KA  © 


]a©k 


nR 


exp  (■ 


-n2R2 


A©K 


A©K 

HI  -  — r-r) 

2n2R2 


(IV.  44) 


This  may  be  simplified  to  the  form 


Qs  =  A2  (A3  +  A4TS°  +  A5q)  (IV. 45) 

A2#  A^  are  constants  throughout  for  a  given  system  with 

fixed  Az  and  A©.  A3  and  Ts°  are  constants  that  are  applicable 
for  each  reactor  element  calculation. 
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Fortunately,  the  evaluation  of  these  constants  is 
simplified  by  the  fact  that  the  series  given  in  equation 
(IV. 44)  converges  very  quickly  and  in  practical  cases  only  the 
first  two  terms  need  be  considered.  It  is  important  to  note 
that  the  constants  need  be  calculated  only  once;  from  then  on 
Q  can  be  evaluated  for  a  series  of  values  of  q  very  simply. 

D 

In  this  later  numerical  confutation,  the  integration  was 
performed  numerically  by  Simpson's  Rule  with  error  correction! 19) . 
The  evaluation  of  the  error  cofunction  is  outlined  in  Appendix  E . 

The  behavior  of  the  gas  phase  is  described  by  a  set 
of  first  order  ordinary  differential  equations,  using  a  time- 
average  value  of  the  solid  surface  temperature  as  the  boundary 
condition.  Thus, 


dTg 

dz 


3(1  -€)  hf 
R  G  C£ 


G, 


N 


m 


Gc 


f  j=i 


dx. 


dz 


AHj 


(IV.  46) 


together  with  the  above  equations  (IV. 3),  .  (IV. 4)  and  the 

appropriate  initial  conditions  in  equation  (IV. 12)  form  an 
independent  system.  The  enthalpy  change  may  be  written  as 


N 


+ 


j=l 


(IV. 47) 
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where  and  c? 

Pi  P2 


are  the  average  molar  heat  capacities  of  gases 


at  Tgi  and  Tg2' 


the  initial  and  final  temperature,  respectively 


IV. 5  Computational  Results 

In  setting  up  a  computational  model  it  is  generally 
desirable  to  proceed  in  a  stepwise  fashion  -  trying  at  every 
stage  to  compare  the  results  either  with  a  known  analytical 
solution,  or  if  that  is  not  possible,  with  a  numerical  solution 
involving  techniques  that  have  been  tested  and  shown  to  be 
valid  and  accurate. 


IV . 5 . 1  Example  Problem  No.  1 

Consider  a  reactor  packed  with  copper  spheres  of 
radius  0.05  ft.  to  be  used  as  a  heat  exchanger  in  an  air 
liquifaction  system.  It  is  desired  to  calculate  the  tempera¬ 
ture  of  the  air  and  of  the  spheres  in  the  exchanger  as  a 
function  of  position  and  time,  under  the  following  operating 
conditions: 


2000  lb/(sq.ft.) (hr) 
4.68  ft. 

280  °R 
350  °R 
0.345 

20  BTTJ/(hr)  (sq.ft.)  (  °R) 

0.23  BTU/(lb)(°R) 

23  lb/(cu.ft.) 


air  flow  rate 

reactor  length 

initial  temp,  of  solids 

entering  air  temp. 

fractional  voids 

convective  heat  transfer 
coefficient 

air  heat  capacity 

entering  air  density 
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This  problem  was  taken  in  part  from  "Applied  Mathe¬ 
matics  in  Chemical  Engineering" ( 20)  and  was  solved  for  four 
different  cases. 

Case  1:  The  results  of  the  analytical  solution  by  equation 
(IV. 19)  for  the  gas  phase,  predicated  upon  a  uniform  particle 
temperature  are  given  in  Table  6a. 

Case  2:  The  results  of  the  implicit  Crank-Nicholson  method, 
outlined  in  Section  IV. 3,  using  a  value  for  K  of  5.4  sq.ft. /hr. 
and  for  ks  of  200  BTU/(hr . f t . ) ( °R)  are  also  shown  in  Table  6a. 

It  is  expected  that  the  two  solutions  would  be  almost  identical, 
since  with  a  high  thermal  conductivity,  (in  fact  copper  was 
used) ,  the  interior  particle  temperature  gradients  would  be 
negligible.  This  can  be  seen  in  Table  6b.  The  close  agree¬ 
ment  confirms  the  validity  of  the  numerical  method  and  the 
accuracy  of  the  finite  difference  expressions  used. 

Case  3:  The  problem  was  redone, using  a  K  of  0.0106  sq.ft. /hr. 
and  ks  of  0.19  BTU/(hr ) ( °R) ( f t ) ,  values  typical  of  refractory 
materials.  These  were  so  chosen  that  significant  temperature 
gradients  in  the  solid  particles  would  exist.  The  results 
of  the  implicit  method  are  tabulated  in  Table  7  for  the  gas 
phase  and  plotted  in  Graph  No.  5  for  the  solid  phase. 

Case  4:  The  decoupling  technique,  by  solving  the  solid  phase 
and  fluid  phase  equations  separately,  was  then  applied  to  the 
problem  presented  in  Case  3.  The  results  of  the  gas  temperature 
are  also  listed  in  Table  7;  but  the  solid  temperature  gradients 
are  plotted  in  Graph  No .  5 .  It  should  be  noted  that  in  the  program 


osotiq  bi  foe  orfi  pnxvloe  yd  pflilqi/ooab  DffT  t f*  esfip  n 

r\ 
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TABLE  6a 

Gas  Temperature  Along  Fixed  Bed  Air 

Liquif action  Heat  Exchanger 
(negligible  intraparticle  temperature  gradients) 


Distance  Gas  Temp  ( °R)  at  Time  0 

from 


Entrance 

(ft) 

Analytical 

Solution 

Implicit 

Method 

AQ=5  sec 

Implicit 

Method 

A0  =  10  sec 

0  = 

100  sec. 

0 

350.00 

350.00 

350.00 

0.267 

337.40 

337  .20 

336.96 

0.609 

321.10 

320.70 

320.31 

1.204 

305.04 

304.43 

304.03 

1.872 

292  .73 

291.95 

291.67 

2  .675 

285.32 

284.58 

284.43 

3.611 

281.82 

281.29 

281.23 

4.680 

280.50 

280.22 

280.20 

0  = 

160  sec  . 

0 

350.00 

350.00 

350.00 

0.267 

341.80 

341.72 

341.57 

0.609 

328.95 

328.78 

328.49 

1.204 

313.63 

313.31 

312.96 

1.872 

299.53 

298.99 

298.69 

2.675 

289.43 

288.71 

288.51 

3.611 

383.73 

283.05 

282.96 

4.680 

281.20 

280.73 

280.70 
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TABLE  6b 

Solid  Temperature  in  Fixed  Bed  Air  Liquifaction 

Heat  Exchange 

(negligible  intraparticle  temperature  gradients) 


Distance  Solid  Temperature  ( °R)  at  160  sec. 

from  by  Implicit  Method 

Entrance 


(ft) 

*2 

t3 

t4 

t5 

0.1337 

328.12 

328.13 

328.14 

328.15 

328.15 

0.4680 

316.74 

316.75 

316.76 

316.77 

316.78 

0.9360 

304.46 

304.46 

304.47 

304.48 

304.49 

1.5377 

293.86 

293.86 

293.87 

293.88 

293.88 

2.2731 

286.49 

286.50 

265.50 

265.50 

265.51 

3.1423 

282.42 

282.42 

284.42 

284.42 

284.42 

4.1451 

280.65 

280.65 

280.65 

280.65 

280.65 
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TABLE  7 

Gas  Temperature  Profile  Along 
Fixed  Bed  Heat  Exchange 

(significant  intraparticle  temperature  gradients) 

th 

At  time  *  5/6  min,  at  the  8  “  increment 


Fractional 

Distance 


Gas  Temp,  ( °R) 


From 

Entrance 

Decoupling 
Method  with 
(E  m  0.01) 

Implicit 

Method 

Simplified 

Analytic 

Method 

0 

350.0 

350.0 

350.0 

1/13 

338,10 

337  .93 

333.71 

2/13 

327.56 

327  .29 

319.80 

3/13 

318.45 

318.11 

308.73 

4/13 

310.76 

310.40 

300.31 

S/13 

304,37 

304.02 

294.12 

6/13 

299.14 

298,83 

289.68 

7/13 

294,92 

294.65 

286.56 

8/13 

291,55 

291,32 

284,40 

9/13 

288,88 

288.70 

282.93 

10/13 

286,78 

286.64 

281.93 

11/13 

285.13 

285.05 

281,26 

12/13 

283,89 

283,82 

280.82 

13/13 

282.45 

282,87 

280.53 

E 


|Qa  -  Qol/ia8  +  Qo‘ 


Where 


Gas  Temperature  Tq  (°R) 
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Graph  I\to 

Tem perature  Gr adients  in  Solids  of 
Regenerative  Bed  Heat  Exchanger. 

(at  1.08  ft  from  inlet) 

I  :  at  time  =  2.0833  min. 
II  :  at  time  =  0.8333  min. 
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A9  =  0.001736  (hr)  (1/48  of  5  min.  cycle) 


so 


A  9  1  /o 

(— ) 1/2 
Ktt 


=  0.234 


and 


R 


pKAQ 


=  11.7025 


hence,  equations  (IV. 43)  and  (IV. 44),  with  the  significant  terms 
remaining  become 


R 


T  = 


2r  J  ttA9K  0° 


r 1 T(r ' ) J  exp 


exp 


m2 


-T 


-(r  +  r') 


4A9K 


dr 


-(r  -  r') 


4A9K 


m2 


RT, 


|  erfc  (yq)  -  erfc  (<5Q)j 


Rq  (  (R  -  r)2 

+  — -  <  ( A9  +  - )  erfc  {yQ) 

r  I  2K 


(R  +  r)  ^ 

( A9  +  - — )  erfc  (6  ) 

2K 


A9  u  2 

-  (R  -  r)  (— )  exp  (-Yq  ) 


A9  i,  2 

+  (R  +  r)  ( — )  2  exp  (-6  ) 

ttK 


(IV. 48) 
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and 


Qs  =  A2  (A3  +  Aa  Te°  +  A,.q) 


4  s 


3Az ( 1  -  £)k 


R 


8.  r  rT(r 1 ) 

0 


KR 


erfc  ( 60)  -  erfc  (y  ) 


where 


6. 


+  p 


-A©  A©  1/2 

-  +  2  (— )WZ 

R  ttK 


+  q 


-A©' 


2R 


4A© 


+ 


AG  1  /9 

(_)  1/2 

irK 


T  \ 


(IV, 49) 


(R  -  r)/2  /a©K 
(R  +  r)/2  /a©K 


For  equation  (IV, 48)  the  temperature  consists  of  three  compon¬ 
ents,  namely,  of  initial  temperature  distribution  T(r) ,  of 
constant  surface  temperature  Ts°  and  of  rate  of  change  of 
temperature  q.  Due  to  the  presence  of  exponential  terms,  the 
numerical  integration  will  not  give  satisfactory  results  unless 
the  radius  is  divided  into  more  than  15  sections.  However, 
the  accuracy  of  this  expression  is  demonstrated  in  Table  8. 

If  T-.  and  are  the  gas  temperatures  at  the 

gl  g2 

beginning  and  ending  of  Az  section  of  the  bed,  then  the  total 
heat  content  change  in  the  gas  phase  for  time  interval  A©  is 


Qg  =  (Tgx  -  V  Ae  6  cf  G 


(IV. 50) 


dr ' 
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TABLE  8 


Temperature  of  Solid  by  Equation  (IV. 48) 


at  T(r')  =  280  ( °R) 

Ts°  =  280  ( °R) 

q  =  0  ( °R/hr) 

as  a  check  on  equation  (IV. 42) 


r 

Temp .  due  to 

T(r) 

Temp .  due  to 

m  O 

S 

Final 

T(r) 

R/20 

279.915 

0.000 

279.92 

R/10 

279.991 

0.000 

279.99 

2R/10 

279.983 

0.000 

279.98 

3R/10 

279.991 

0.000 

279.99 

4R/10 

280.008 

0.000 

280.01 

5R/10 

279.971 

0.020 

279.99 

6R/10 

279.572 

0.436 

280.01 

7  R/10 

274.772 

5.220 

279.99 

8R/10 

245.744 

34.278 

280.02 

9R/10 

153.049 

126.924 

279.97 

R 

-14 

1.3  x  10  ^ 

280.00 

280.00 
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The  calculation  procedure  may  be  stated  as  follows: 

1.  estimate  the  surface  temperature  from  equation  (IV  .6)  so 


12Ar  hf 

(  - )  t  +  18T.  -  9T,  ,  +  2T 

9l  n 


k. 


■h-1 


n-2 


(11 


+ 


12  Ar  hf 
- - ?_) 


(IV.  51) 


where  Tn  is  the  outermost  solid  temperature* 


2 .  obtain  q  as  a  first  approximation 

m  m  o 


q  = 


A© 


3.  calculate  Q  from  equation  (IV. 49) 

4.  solve  equation  (TV. 46)  with  reaction  term  neglected  to 

9et  Tg2 

5.  calculate  Qg  from  equation  (IV. 50) 

6.  if  |  Qs  -  Q  |  ^  A,  a  given  allowable  difference  proceed  to 
another  increment  and  repeat  the  whole  process  of  computation, 

7.  if  |qs  -  Qg|>  find  a  new  <3  from 


q 


(Qs  +  Qg) 

2 


a2  a5 


(IV. 52) 


and  go  to  3 . 
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It  can  be  seen  from  Table  7  that  the  close  agreement 
with  the  implicit  method  suggests  the  calculational  techniques, 
employed  in  the  decoupling  method,  are  free  of  gross  errors 
and  serve  as  a  check  on  the  soundness  of  the  mathematical 
model . 


IV  . 5 . 2  Example  Problem  No.  2 

The  decoupling  method  has  been  proved  to  be  prac¬ 
ticable  and  accurate,  so  that  it  may  be  applied  to  the  design 
of  a  regenerative  fixed-bed  reactor.  The  problem,  with  design 
data  taken  from  Kumagi  et  al(21),  is  stated  as  follows: 

A  regenerative  refractory  furnace  is  used  to  produce 
acetylene  by  the  thermal  cracking  of  methane.  Pure  methane 
at  950°K  is  fed  at  a  rate  of  20  lb-moles/( sq.ft . ) (hr)  into  a 
pebble  bed  of  alumina  with  radius  0.5  inches  and  a  uniform 
initial  temperature  of  2000 °K. 

The  main  reaction  is 


2CH4 


c2K2  +  3H2 


(IV. 53) 


where 


AH 


1400 °C 


exp  (12.836  -  19214/T  )  sec'1 

y 

exp  (4.358  -  5992/Tg)  sec-1  (Tg  °K) 

26,654  BTU/lb-mole  CH4  converted 


The  bed  is  3  feet  in  length  and  has  a  porosity  of 
0.4.  The  convective  heat  transfer  coefficient  between  the  gas 
and  the  solid  particles  is  estimated  to  be  46.1  BTU/(hr)  (sq.ft .)(  °R)  . 
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Other  physical  data  are  obtained  from  Perry's  "Chemical 
Engineers'  Handbook". 

It  is  desired  to  find  the  temperature,  fractional 
conversion  and  the  pressure  drop  at  various  locations  along 
the  reactor  as  a  function  of  time. 

Accordingly  the  computational  procedures  follow 
those  outlined  in  Case  4  of  the  previous  problem,  with  the 
inclusion  of  reaction  and  pressure  drop  in  the  gas  phase. 

Since  the  reaction  is  given  to  be  reversible  and  1st  order, 
if  ideal  gas  behavior  is  assumed,  then  the  rate  equation 
can  be  written  as 

P  ,  'I  lb-mole 

R(x,P,T  )  =  -  Ik^l  ~  x)  ~  k2X//2 - (IV. 54) 

R*T  (1+x)  '  sec. ft2 

9 


The  fractional  conversion  becomes 
dx  P  3600 

dz  Gm°R,T  (1+x) 


exp 


(12.836  -  19.214/T  ) 


—  •  exp 
2 


(4.358  -  5.992/T  ) 

SJ 


(IV.  55) 


Similar  to  that  for  moving-bed  reactor  the  energy 


balance  in  gas  phase  is 


dT 


dz 


2  = 


kf  as 


G, 


(Ts  +  qA9/2  -  Tg)  - 


CP2  3CP3 

( -C  +  - -  +  - -) 

P1  2  2 


m 


— 

dx) 

x  3x 

(T  -  1673)  +  26,654 

— 

/ 

(1  -  x)  C  +  —  C  +  —  C 

g 

dz 

P1  2  p2  2  p3 

(IV.  56) 


.  "  jCc  odhnsU  *8’  o 
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At  1000 °C  and  35%  conversion,  which  represents  the 
average  values  of  the  parameters,  the  viscosity  of  the  gas 
is  0.03  cp .  The  density  is  given  by 

16 

pg  =  - 

(1  +  x)  (359)  (T  /273)  (14.7/P) 


0.83  P 


(1  +  x)  T 

g 


(IV.  57) 


Upon  substitution  into  equation  (IV. 4),  the  pressure  drop 
becomes 

dP  (1  +  x)  T 

—  = - ^  x  8.35  -  10"4  (psi/ft)  (IV. 58) 

dz  P 


The  results  are  all  plotted  in  Graphs  No.  6  and  7. 
There  is  no  experimental  nor  calculational  check  on 
these  results,  although  the  general  trend  appears  reasonable. 


IV. 6  Summary 

1.  For  fixed  bed  heat  exchangers  with  no  chemical  reaction, 
the  proposed  implicit  method  which  involves  dividing 
the  solid  particles  into  a  number  of  incremental  shells 
and  using  a  time  average  for  the  solid  temperature  and 
a  position  average  for  the  gas  temperature  is  a  general 
and  stable  method  that  should  be  satisfactory. 


Gas  Temperature 


Graph  No.  6 
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Gas  Temperature  Profile  along  the  Fixed  Bed  Reactor 
for  the  Production  of  Acetylene. 


Length  of  Reactor  z  (ft)  . 


Fractional  Conversion 


Graph  No.  7 
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The  Fractional  Conversion  of  Methane  along  the 

Fixed-bed  Regenerative  Reactor. 


I  :  at  time  =  0.5  min. 
II  :  at  time  =  3.0  min. 


Length  of  Reactor  z  (ft) . 
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2 .  For  fixed  bed  regenerative  chemical  reactors  with  reaction 
in  the  fluid  phase  only,  the  very  large  problem  involving 
a  system  of  non-linear  partial  differential  equations  can 
be  reduced,  by  a  decoupling  technique  into  two  much  smaller 

problems  involving,  first,  a  linear  partial  differential 
equation  and,  second,  a  set  of  first  order , non-linear , 
ordinary  differential  equations.  An  efficient  computational 
technique  for  solving  the  former  is  presented;  conventional 
explicit  methods  are  employed  on  the  latter. 


1-71 


NOMENCLATURE 


G 

m 

Gs 


AHTfa 

k 

ks 

K 

P 

QG 

Qs 

r 

R 

R' 

R(x,P,T  ) 
y 


Total  surface  area  per  unit  volume  of  reactor  space 

Concentration  of  component  a 

Molal  heat  capacity  of  gas 

Molal  heat  capacity  of  component  i 

Specific  heat  of  gas 

Specific  heat  of  solid 

Diameter  of  solid  spheres 

Inner  diameter  of  reactor 

Mass  flow  rate  of  gas  per  unit  area  cross-section 

Molal  flow  rate  of  gas  per  unit  cross-section  area 

Mass  flow  rate  of  solid  spheres  per  unit  cross-section 
area 

Heat  transfer  coefficient  between  gas  and  solid  spheres 

Heat  transfer  coefficient  between  gas  and  reactor  wall 

Heat  of  reaction  at  reference  temperature  T^ 

(positive  for  endothermic  reaction) 

Reaction  rate  constant 

Thermal  conductivity  of  solid  spheres 

Thermal  diffusivity  of  solid  spheres 

Total  pressure  at  any  point  z  in  reactor 

Enthalpy  of  gas 

Heat  content  of  solid  spheres 

Radius  variable  from  center  of  solid  sphere 

Radius  of  solid  sphere 

Gas  constant 

Rate  of  conversion  of  limiting  reactant  per  unit 
volume  of  reactor 
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w 


x 


Time  variable 

Temperature  variable  of  solid  sphere 

Temperature  variable  of  gas 

Surface  temperature  of  solid  sphere 

Temperature  of  reactor  wall 

Velocity  of  solid  spheres  in  moving  bed 

Fractional  conversion  of  limiting  reactant  at 
point  z  in  reactor 

Distance  variable  along  reactor  from  the  entrance 

Fractional  void  in  bed 

Density  of  gas  at  £>oint  z  in  reactor 

Density  of  solid  spheres 

Reduced  time  variable  defined  in  equation  (IV. 9) 
Viscosity  of  gas 


Superscript 


initial  condition 
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APPENDIX  A 

Finite  Difference  Expressions 
Which  are  Third  Order  Correct 


If  Tn_2'  Tn-1'  Tn'  Tn+1  are  known  the  spacings  a ,  b,  c 

are  not  necessarily  equal,  then  by  Taylor's  expansion 


?  , ,  -  T 

n+1  n 


a  T'  + 
n 


a 

2 


ip  i  i 

n 


___  ip  i  i  i 
6  n 


T  1 

— 

T 

— 

-b 

T' 

+ 

— 

«P  i  i 

_ 

— 

«P  i  i  i 

n-1 

n 

n 

2 

n 

6 

n 

c2 

C3 

T  0 

— 

T 

-c 

T' 

+ 

— 

m  «  • 

■  — 

rp  1  I  1 

n-2 

n 

n 

2 

n 

6 

n 

Hence, 


T'  = 
n 


2,  2 


2  2 


2,  2 


(Tn+i  -  Tn) c  b  (b-c)  +  (Tn_1-Tn) a  c  (c+a) - (Tn_2-Tn) a  b  (a+b) 


T  •  i  — 

n 


(Tn+1-Tn)bc(b2-c2)+(Tn 1-Tn)  ac ( a2-c2)  +(Tn 2-Tn)  ab(b2-a2) 

A/2 


where 


.  ,  ,  2  2.  2  2.  ,  2  ,2 
A  =  abc  J  a(b  -  c  )  -  b(c  -  a  )  -  c(a  -  b  ) 
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Differential  Finite  Difference  Expression  Geometrical 
Form  Representation 


2T  -,  +  3T  -  6T  ,  +  T  0 
n+1  n  n-1  n-2 

T  O 

n-2 

Tn-i 

T  T  , 

n  n+1 

£ • 

dr 

6  Ar 

..2 

a  Tn 

Tn+1  "  2Tn  +  Tn-1 

Sr2 

Ar2 

32Tn+%  -  15Tn  -  20Tn_1+3Tn_2 

Tn-2 

Tn-i 

T  T 
n  n+k 

dr 

30  Ar 

d2T 

n 

16Tn+^25Tn+10Tn-l-Tn-2 

dr  ^ 

5  Ar2 

11T  -  18T  ,  +  9T  ,  -  2T_  , 

n  n-1  n-2  n-J 

Tn-3 

T_  9 

n-2 

TV  TV 

n-1  n 

dr 

12  Ar 

-.2 

a  Tn 

2Tn  -  5Tn-l  +  4Tn-2  “  Tn-3 

dr2 

Ar2 

9Tn+^ 

184Tn+Js  -  225Tn+50Tn-l-9Tn-2 

Tn-2 

T  , 
n-1 

T  T  ,  j 
n  n+% 

- - - * - 

dr 


60  Ar 
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APPENDIX  D 

Outline  of  Runge-Kutta  Method 


dT 


dz 


=  (z;  T±t  ...  ,  T±t  ...  ) 


k. 

1 


=  Az  T,  (z;  ...  f  Ti#  ...) 


k.  '  ' 

1 


Az  k-|  k  f 

Az  T,  (z  +  — ;  T1  +  — ,  ...  ,  Ti  +  •  •  •  ) 


k.'  ‘ 
i 


Az  k ' •  k! • 

,  1  i 

Az  T.  (z  +  — ;  T  -h  ,  ...  ,  T.  +  — ,  ...  ) 
1  Z  \  Z  1  z 


klv 

i 


Az  T .  ( z  +  Az ;  T1  +  k  '  '  ' ,  ...  ,  T  -I-  k '  '  ' ,  ...  ) 

1  1  1  i  i 


Therefore, 


AT. 

i 


1 

6 


k!  +  2  (k !  ' 

i  i 


-i-  k !  '  ' )  +  k 

i 


n 


All  the  k|  are  evaluated  before  the  kj '  are  cal¬ 
culated,  similarly  the  k'1'  and  k1V  are  evaluated  sequentially 

i  i 

If  Runge-Kutta-Gill  is  used  the  coefficients  will  be  dif¬ 
ferent  but  they  are  given  in  reference  (8) . 
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APPENDIX  C 

Stability  of  the  Diffusion  Equation 


Consider  an  I.V.P.  of  heat  conduction 


with 


du  d^u 

—  =  K  -  0  <  x  <  L 

cO  dx^ 


u(x,  0)  =  ^(x) 

u  ( 0 , 9)  =  0 

u(L, 9)  =  0 


(1) 


(2) 


If  u(x,9)  is  the  exact  solution  and  u-  is  the 

J  / n 

finite  difference  solution  then 


u(x,9) 


I 


nirx  nir  2 

A  sin  -  exp  (-K  ( — )  9) 

L  L 


(3) 


where  AR  is  obtained  from  expanding  V (x)  into  an  infinite 
orthogonal  series  such  that 

L 

2  p  m t 

A  =  —  /  (x)  •  sin  (  —  •  x )  dx 

n  L  J  L 

o 

Expressing  (1)  in  finite  difference  forms  of 
(i)  explicit,  forward  difference 


u 


j  /  n+1 


=  u 


‘D#n 


P 


2u 


u  .  . ,  ) 

1+1/ n' 


(4) 
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(ii)  implicit,  backward  difference 


“uj-l,n+l  +  ^  u  j , n+1  u j+1, n+1 


=  u 


D#n 


£ 


(iii)  implicit,  Crank-Nicholson  method 


■u  .  ,  ,,  +  2(1  +  — )  u .  . 

3+1, n+1  p'  3, n+1 


u j+1, n+1 


u  .  -i  _  -  2(1  -  — )  u .  +  u.., 

D-l/H  p'  3 , n  J+l,n 


(6) 


where 


(3  =  AG/Ax' 


Now  define  the  error  as 


€  j 7 n  =  uj/n  "  u(nAx'  3A9) 


and  substituting  into  (4)  it  gives 


€.  =  P  €  .  _  +  (1  -  2p)  €  .  +  p  €  .  _  (7) 

3, n+1  3+lfH  3  # n  3-1, n 


By  separation  of  variable, solution  of  equation  (7)  may  be 
written  as 


-  _  ,  a3A9  iYnAxx 

€  •  =  f  (e  J  ,  e  ) 

3 , n 


(5) 


(8) 


r 
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Substituting  into  (7)  and  simplifying  gives 


aA9 

e 


P  e1^  +  (1  -  2P) 

1  -  4f3  sir/  (-y) 


+ 


pe-i^Ax 


(9) 


In  order  that  the  error  € .  at  any 

jtn  * 

will  not  grow  as  ©  increases,  it  is  necessary 
that 


fixed  point,  nAx, 
and  sufficient 


or 


(10) 


By  this  procedure  one  can  establish  the 
ratio  of  the  above  finite  difference  forms: 

(i)  explicit 

a  A© 
e 


2  y  Ax 

^1-4(3  sin  ( - )  ^  1 


stability 


(ID 


(ii)  backward  difference 

a  A© 


9  yAx  . 
^  1/(1  +  sin2  - )  ^ 


(12) 


(iii)  Crank-Nicholson  form 


4f3 


.  2 
sm 


yAx 

( - ) 


2 


2  ^Ax  / 

/(I  +  4p  sin - )  1 

2 


(13) 


Thus,  (ii)  and  (iii)  are  unconditionally  stable. 
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APPENDIX  D 

Verification  of  Equation  (IV, 43)  in 

Satisfying  the  Initial  Condition 


Since  Equation  (IV. 43)  is  the  sum  of  equations 
(IV. 31)  and  (IV.  42),  at  9=0  equation  (IV. 31)  is  zero;  it  is 
necessary,  therefore  to  prove  that  equation  (IV. 42)  does 
satisfy  the  initial  condition  T(r) . 

Lemma: 

If  f(x)  is  an  even  function  of  x,  which  can  be 
expanded,  so  will  be  f(x"t.2nR),  in  a  Fourier  series  of  cosines 
of  multiples  of  ttx/R,  then 


oo 


oo 


)  f(x+2nR)  =  — 

L>  R  u 

n=-oo  0 


f(x)  dx 


2 

+  — 
R 


oo 


n=l 


oo 


cos 


mrx 


R 


mrx ' 

f(x')  cos  -  dx' 

R 


(i) 


0 


provided  the  integrals  are  convergent  and  the  series  will 
converge . 

Proof  of  the  lemma  can  be  referred  in  any  standard 
text  on  Fourier  Transforms. 

-x2/4KO 
e 


Now  let  f(x) 


,  then 


_ 
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GO 


n--oo 


exp  { 


(x  +  2nR) 2 


4K0 


1 

R 


co  2 

n  X 


0 


exp  (-  - )  dx 

4K0 


R 


CO 


n=l 


COS 


mrx 


R 


co 


0 


x“  mrx 1 

exp  ( - - ■)  cos  — - dx* 


4K9 


R 


JttKQ 

R 


1  -!-  2 


CO 


t  J 

n=l 


mrx 

cos  exp  ( 

R 


2  2 

Kn  TT  0 


)  f  (ii) 


R" 


Here  letting  fj  =  x/2  KQ  and  b  =  mr/KQ/R 


0 


co  p 

p  -X  mix 

exp  ( — — )  cos  - -  dx  = 

4K0  R 


co 


2  J KQ  j  e 
0 


■r 


cos  (2brj)dr|  (iii) 


From  the  contour  integral 


z‘ 


dz 


0 


for  z 


1  + 


(iv) 


where  c  is  the  contour  of  the  rectangle  -a  <C  IJ  <  a,  0  ^  ^  ^  b 
and  taking  the  limit  of  a  approaches  infinity  it  yields 


co 


-»7 

e  *  cos 


( 2b  »|  )  dcj 


1/2  •  /  TT  exp  (  -d  ) 


(v) 


0 


Hence,  the  equation  (ii) 


1-83 


Replacing  x  by  (r1  -  r)  and  (r  +  r)  ,  equation  (IV. 42) 
yields  upon  simplification 


u 


„  R  oo  2  2 

2  p  q-1  Kn  7T  9 

—  /  T(r')r‘  )  exp  ( - - - ) 

rR  J 

0  n=l 


R“ 


nTrr  nTr 1 

sin - sin  -  dr' 

R  R 


(vi) 


which  is  equation  (IV. 32). 


at  9  =  0 


u, 


1 

R 


oo 


n=l 


sm 


mrr 

R 


R 


r/2 


0 


mrr ' 

r'T(r')  sin  -  dr 

R 


(vii) 


Mow  assuming  that  rT(r)  can  be  expressed  into  an 

mrr 


infinite  orthogonal  series  of  sin 


R 


l  .e . 


rT(r)  = 


oo 


n=l 


Bn  sin 


mrr 


R 


then 


R 


0 


r,rr(r')  sin(^—  r  1 )  •  dr  1 


n 


R 


o  mr 

sin^(  —  r  1 )  dr  ' 
R 


0 


' 
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R 


r'T(r')  sin( 


mr 

R 


)  dr 


R/2 


Therefore, 


rT(r) 


nir 

sin( —  r) 
R 


R 


R/2  d 


r 1 T(r ' ) 


mr 
sin(  — 


R 


so 


u  =  T(r) 
o 


Q.E.D. 


(viii) 


')  dr '  (ix) 
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APPENDIX  E 


Numerical  Evaluation  of  Error  Cofunction 


By  definition 


and 


erfc(x)  =  1  erf(x) 


erf(x)  sb 


0  x  ? 

f  e ”t  dt 


ft 


0 


From  mathematical  tables  values  of  erf(x^)  are 
read  into  the  memory  of  the  digital  computer.  To  evaluate 
erfc(xj)  where  x^ <  Xj  <  x^+]  ,  let 


Ax 


X  .  -  X  • 

_ i 


and 


^k 


exp  (-(x^  +  kAx)  ) 


for  k  =  0,1,  . . . ,  6 


Then  by  Newton's  forward  formula  of  integration  for  the  6th 


degree 


1  -  erf(x^)  - 


2  Ax 


140  ]~¥ 


(41y  +  216y  +  27yn  +  272y 

o  1  ^  j 

+  27y .  +  216y  +  41y,) 


orfc(Xj) 


1-86 


APPENDIX  F 

Listings  of  Computer  Programs  and 

Sample  Outputs 

1.  Moving  Bed  Heat  Exchanger  Design 

2 .  Moving  Bed  Reactor  Design 

3.  Regenerative  Bed  Heat  Exchanger  -  Analytical 

4.  Regenerative  Bed  Heat  Exchanger  -  Implicit 

5.  Regenerative  Bed  Heat  Exchanger  -  Decoupling 

6.  Acetylene  Production  in  Regenerative  Bed  -  Decoupling 


5  JOB 
ST  IMF 

S  I  B JOB  MOVING 
SIBFTC  HFAT 
C 


1  5  00  0  3  MOVING  F  r  D  HEAT  "  X  C !- !  A  N  C  C  R  • 

5  *  5  0  0 
GO 

NODFCKt NOL 1ST 


r 

r 

r 

r 

r 

r 

c 

c 

c 

c 

c 

c 

r 

c 

c 

r 

c 

r 

C 

C 

C 

C 

c 

c 

c 

r 

c 


i 

? 


OFSIGN  OF  MOVING  BFD  H^A  T  EXCHANGE0  BY  FXPRESSlNG  THE  DIFFUSION 
FOUAT I  ON  I N 1 0  A  T  '  .  .  •  .  ,  '  ,  h  1 

GA  S  FOLIATION  ARF  SOLVFD  r  R(  NGF-KUTTA-GILl  HOO* 

All  INPUT  hai  a  A  R  c  IN  LB,  mq,  py,  BUT*  DFG  F*  ! 

M A  ••  i  10  •  OF  DIVISIONS  IN  RADIUS  SOLID  i  •  • .. 

IF  NFCFSSARY,  CHANGF  T  HF-  FORMAT  FI  PUT  . 

*  fi  ■  f  solid  and  gas  r f s p f c t i v f l t • 

CS*  CP  -  I  I  SOLID  AND  •  '  I  .  l  . 

:i  ’  IGO  -INI  L  TEMPERATURES  OF  SOI  ID  AND  3A5 
nS  -  OFNSl  sOLID  SPHER I  • 

R  -  RADIUS  OF  SOLID  SPHERE  IN  INCHES. 

SK  -  THFPMAL  CONDUCTIVITY  OF  SOLID  SPHERES. 

VD  -  FRATTIONAL.  V^fi)  OF  THE  MOVING  EFT). 

HF  -  HFAT  TRANSFER  rOPFFiriFNT  BFTWbRN  GAS  ANn  rOLin  ,'PH;rPFc 
Dt  -  INCRFMENT  OF  FXCHAN( 

M  -  NO.  OF  DIVISIONS  1  •  '  :  i  '  . " 

TS  -  TEMPERATURE  at  THE  SURFAC 

RTS  -  NORMALIZED  SURFACE  TFMPERATURE  ( TS-TO ) / ( TGO-TO ) . 

IG  -  GAS  TFMPFRATURF  VARIARLF. 

PTG  -  NORMALIZED  GAS  TEMPERATURE  ( TG-TO )/( TGO-TO) . 

OUTPUT  DATA  ARE  IN  SETS  OF  TH9FF  LINFC 

DISTANCE  IN  FEET  AND  REDUCED  DISTANCE  OF  TH*D I  ST / ( US*R*R ) 
TS,  RTS*  TG  AND  RIG 

SOLID  TFMPFRATURFS  AT  ALTERNATING  SPHFRICAL  SHELL  FROM 

DTMFNSTON  F(?0),  Q  A  (  2  0  )  »  QB ( 2  0 )  *  OC(20),  00(20),  (20),  Y(20) 

PFAD  (5,1  )  GS,  CS,  DS,  SK ,  R,  VO 

FORMAT  ( F  F .  2  ,  FF.b,  F 7.2,  F7.3,  F5.2,  F5.2) 

RFAD  (5,2)  GA,  TP,  HF,  DL,  M,  TO,  TGO 

FORMAT  (F8.2*  F  6  •  3 ,  F7.2,  F6.3,  2X,  12,  E7.2,  F7 • 2 ) 

M A  =  M+ | 

M0=M- ] 

F  m  =  M 

R  =  R / 1 2.0 
H=  R/FM 

T H  =  S K  /  (  CS*l)S  ) 
u s  =  g s / ( ns* ( i .  i-vn ) ) 
ap=Tm/( m*h*US ) 
p  ED=  TH/ ( US*R*R ) 

)A=^,iJ*H|-*(  "I  .  0  —  V I )  )  /  (  o  A  *  C  P  *  i '  ) 

S] =12.U*FM-1 .8.0 
S  2  = 1 2  0 . 0  *  F 5  —140.0 


3 . 

• 

•  • 

•  • 

• 

• 

•  ’ 

0 

0 

. 

1 

• 

V-, 

. 

<  • 

* 

• 

• 

. 

• 

• 

—  1 

• 

• 

• 

• 

• 

«  • 


•  t 


( 


.  ( 

0 

(  r , 

•  <  . 

'AM '  ’  n  1 

0 

( <;* 

I  — 

»  •  _  1  <  3 

.S  r\  V  =  « 

M  ^  V  >^H 

(  > )  \  *  -'  •  r 

(  i  1  —  .  r  )  * 

(  -»  |*M'fM)\UT3<lA 

.  r  )  >  nf 

.  -  .  1  ~  ’ 
h  *  .  v r  =  <; 


53  =  3  30 . 0*FM-90. 0 

54  = 1 9 2 . G*  F M+ 32 . 0 

S  5  =  3  0  .  U  *  (  2 . 0  *  h  M  -  1  .  0  ) 

X  A  =  SORT ( 2 . U  ) 

P  A  = ( X A- 1 • 0 ) / 2  *  0 
PR=  (  2  • 0-X A ) / 2  •  0 
Pf=-XA/2.) 
pn=  (  2  .  )+X  A  )  /  2  •  ) 

RT=AU.0*HF*H/S< 

R I T  =  1 P4.0  +  R  r 
dtst=o.oo 

no  2  T  =  1  ,M 
3  T (  I  ) =  TO 

T  (  M  A ) =  T  GO 

1 21  DO  10  r  =  1  » MA 
iu  y ( i )  =  t ( n 
C  R  =  1  .  0 

100  Ml  )  =  AP*  (  -  ?6. 0*Y  (  1  )+27.0*Y(2)-Y(3))/R.0 
F(?)=AP*(-Y(1 ) -1 2 . 0*Y ( 2 ) +1 3 . G*Y ( 3 ) ) /9. 0 
no  ii  J=R ,MP 
F  j  =  j 

F  A  =  2  •  0  *  Y  (  J  —  2  )  /  (  A  •  0 * F  J”  3  •  0  )  +  (  2  •  0* h  J—  ^  •  0  )  •#  y  (  J  —  l  )  /  (  2  •  0*  F  J— " 
B  =  -4*0*( FJ-1 .0 )*Y ( J  )  / ( 2.  J-l • 0 )  +  ( 6 •  •  . 

11  F ( J ) =AP*( FA+FB ) 

F  S  =  (  B  T*  Y  (  i)+  .  v  )  - ;  .  ( M- 1  )  ■»  ) .  0*  Y  ( M-2  )  )  /  B  T  7 

F ( M ) =AP* ( -Si  * Y ( M-2 ) +S2*Y ( M-i  ) -S3*Y ( M ) +  S4*T S ) /S5 
F ( MA ) =  BA* (  Y ( MA ) -IS ) 
l  i-  (  •  G  T  •  .  U  )  GO  i  0  !  3 

DO  12  I =1 *MA 
Q  A  (  I  )  =  F  (  T  )  *01 

12  Y (  T  )  =  T (  T  ) +QA (  T  ) / 2 . J 
fp=CP+l .0 

GO  TO  1 00 

13  If-  (OP.GT.2.G)  GO  TO  IF 

DO  14  1=1 *MA 

OR (  1  ) =F (  1  ) *OL 

14  Y  (  I  ) sT (  T ) +QA ( 1  ) *PA+QB i I ) *PB 
CR=( R+l  .0 

GO  TO  100 

15  IF  (CP.GT.3.0)  GO  TO  17 

no  1  A  1  =  1  ,MA 
OC  (  I  )  =  F  (  I  )  *D|. 

1  A  Y  (  T  )  =  T (  1  ) +QR (  1  )  *pr  +QC {  1  )  *PO 
rp  =  rp+i  . 0 
GO  TO  100 

17  no  lfl  r  =  1 ,ma 

on (  I  )  =  F ( I  ) *OL 

13  T  (  I  )  =  T  (  I  )  +  (  QA  (  1  )  +Qi)  (  I  )  )  /  6 . 0+  (  p  '*08  (  I  )  +PD*0C  (  I  )  )  /  3  •  0 
F S= ( B T *1  ( MA ) +2 2  5 • 0* T ( M ) - 5 0 • 0* 1  (  • - '  ) +9 •  (  - / )  )  / 
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(  I  )  '*  (  I  )  Y 

.  r  =  °  ■> 

.  i  i  ■  ,m  +  i  )  .  $  -  (  r  )  -  )  *cja  -  f  c  ) 

,  F  -( 

L  =  LO 

.  .  .  .  .  .  • 

i  )  *  <4  A  r  (  t  )  4 
+•  (  ,  M  )  T  *  R  )  =  <L  ' 

R  a  (  A  *•'  ) 

(  •  I  •  ^  »  <30  ) 


in*<  j)  R*  ( T  )  am 
.  c  \  (  t  )  '(>  +  (  [  )  r  -  f  T  )  Y  c 
.  r  +  n  -  o  -> 

(  \'0) 

/■  .  -  f  4-N r  on 


J(i*  <  I  )  -*=  (  I  )  RO 

>+  «-.*  <  r  )  >  (  r  J  =  (  I  )  r 

(  •  c  .  >  •  o  > ) 

•  ■  .  '  =  j 

m*  (  t  )  Ris  (  t  >  ")C> 

- 

. r +oirin 
•  f  •-  I  r  : 


]  -  MS 


R TS= ( TS-T 0 ) / ( T bO-T  0 ) 

R  T  b= (  F ( MA ) -TO ) / ( T GO- T 0 ) 

I)  I  ST  =  I>  I  ST  +  OL 

P  F  I)  I  57  =PF  !)*i>  I  ST 

WR  T  TT  (  F  9  5  )  O  T  8  T  -*  PFD  T  ST 

FORMA  i  (  1  Flu  9  -JXf  11  HOT  ST  ANTE  =  *  FF.8*  4H  FT.  9  5X  *  i4H:;  )  >  , 

^  '  H  »  ?X  !  F  1  8  .  F  ) 

W  R  T  T  F  (  F  *  F  )  i  »  R  I  9  •  (  m  A  )  9  R  T  G 

h  FORMAT  <1H  9  ?X*  5HTS  =  9  •  2  9  4H  OR  *  r  1 r  •  F  9  "  X  ♦  FHTG  =  9  *-  7 .  ;  , 

] 4H  OR  ♦  h  1 5 . F ) 
w  R  I  T  F  (6*7)  ( T (  I  ) 9  I  =  1  *  M  9  2 ) 

FORMAT  (IN  9  F  x  9  FFil.2) 

IF  (9TG.LT. 2. 5)  GO  TO  101 
WRITE  (698) 

8  FORMAT  (IF)- 9  18X9  2QHF'MO  OF  CALCULATION.  ) 


FND 


$F/\|TRY 

1 2  8  .  -  .  2  0  C 

2500.00  0  .  ?  c  0 


hfa  t 

1  F2  .  c  1  .  ^uo  1.00  0.8' 
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15uu03  MOVING  B  i- » > 

I )  I  S 1  ■  =  0  .  1  0  0 

T  S  =  1  5  3  .  0  ] 

TOO. 00 

O I  S  T  A  N  r  F  =  0  .  ?  C  0 

TS  =  1  7?.  1  7 

1  u  u  •  1 0 

Ml  •  -  =  0.3 

TS  =  167.43 

1 u  u • Q6 

I)  I  S T A N C F  =  0.400 

TS  =  201.04  OP 

103.48 

=  0.5 

IS  =  210.77 

] u  8 . 02 

OISTANCF  =  0.600 

TS  =  225.96 

114.39 

0 1  STAN  = 

237.83 
1  2  2  .  J  9 


D  I  S  T  A  N  (  u  =  0 . 8  C  0 

IS  =  249. 48  OP 

131.04 

i )  i  S  FAN  =  i-  •  ■  >  0 

l  a  260. ■ 

1  4  u  •  b  4 

0  1  S  i ANC  E  =  1.00  0 

T S  =  272.3 

150.75 

D I STANCF  =  1.1 00 

7s  -  283.73 

161 .22 

o I  =  1.260 

I  s  =  -1  u  1.79 

178.43 


( SAMPLE  OU FpUI ) 


I  N . 

FT.  REDUCFD  LENGTH 

•  ]  4  7  F + 0  C  T  G  — 

10'  .01  1 00.1  - 

f  1  .  REDUCED  LFNGl  H 

0 . 7  2  1  7  3  7  F  4-  0  0  Trn  = 

1  0  0  .  3  9  I  0  2  •  4  1 

'  .  1  )  U  C  E  0  LENGTH 

0.874  31  8  E  "P  0  0 

102.19  107.12 

F  1  •  RE  OUT  El )  LE NG  I  H 

J.1U1041F  +  01  7  G  = 

1  u  6  .  u  5  114*21 

FT*  REDUCED  lfngth 

0 .  I  1  8  7  6  7  0  +  01 
1  1  1.89  1 2  2 • 7  3 

r  7  .  C  E  0  L  E  N  G  T  H 

0 • 1 4 5964 F +01 
119.31  1 3  2 • i 6 

■  ; .  REDUCED  LENG CH 

0.1  x7629E  +  0l  7  G  = 

1 2  7 • 8  9  1  42.1 


.  PFDU<  Lf 

0 . 1  4  9  4  7  7  F  +  0  1  1  G  = 

1  3  7.30  1  5 2  •  6  0 

r  I  •  REOUC t1  D  LENG T rt 

'j  •  160980E  +  01 
147.29  l63.2o 


r-  i  .  REDUCED  LENGTH 

. i 72 3 8 7E+(  t  = 

15  7.67  1  7  4 . l 3 

F  .  .  PFDUCFD  L  F NG T H 

0.1 83728E+0 
1  6  8  •  7  2  185.10 

FT.  EDUCED  LENGTH 

0 • 2  C 1 7  6  6  F  +  0  ] 

1  8  .  7  1  2  02.50 


0  •  1  4  7  9  o  4  r  —  0  1 

215.  0.11:922E+01 

1  0  2  •  2  -  i  1  7  .  ~  ~ 

0 . 2  °  9  m  0  s  c  _  -  ] 

2  -  •  •  0.12 
1 10. m  1  136.17 

0 • 44 9 7 1 2  F-o  1 

241 .  1  0 • 141606E+  G  1 

1  2  1 . 2  0  1  5  1  •  "  1 

0.59961 6  F  —  0 1 

2  •  3 • 8  0 •  1  5 3627F+Q 

132.2  163.63 

0  •  7  4  9  "  2  0  F — o  1 

2  5.  . 

1  4  7  •  7  7  1  '  -  •  6  0 

0 •  -  7  942 3  F  —  0  1 

277.  0 . 1 7  7 1 2  5  E+ C 1 

1 3 :  • 3  7  190.78 


0.1  ’)  4  9  3  3  F  —  7  0 
28  8.  O.i  86  56"^E+C 

1  6  6  •  /  2  0  2  .  J  r 


0.1  1  9  9  2  3  F - 0 0 
-  •  0 . 1 9  9  9  3  0  F  +  C 

176.  1  2 1 4 •  7  0 

0  •  i  3  4  9  i  4  F  -  0  0 
311.24  0.21 1241E+( 

169.41  226.26 

0. j  49904 F— 00 

-  2 .  .222  4 

20'J.6»  2  3  7.^1 


0. 1 64P94F-00 

. 

211  .  73  2  4  9 • 0  8 

•  7  9  E— C 

351.7  .251” 

2  2  7.91  26  7  .  1  7 
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SJ06  15u001  MOVING  6£d  RPaCIOR  DESIGN* 

5b  T  I  ME  5  »  5  0  0 


$  I  B  JOB  MOVING  GO 

5b  I  ^  r  T  (  Rt-ACl  NuL  I  S  T  *  NODF  c~< 
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A  SPFGJF1C  PROGRAM  DFVFLOPED  IN  V I FW  r » -  ThF  MOVING  PFD  HP  AT 

pxrHANGEP  for  Thermal  fracking  of  n-putanf. 

All.  T  N  PUl  OATA  AP(-  in  I  R  *  HD*  FT*  HTU*  DFG  F*  UMl  ESS  SPR  C  T  h  I F  f>  # 

the  radius  of  solid  Sphere  has  bfen  divided  into  5  fqual  segment*. 

dZ  —  INCREMENT  OF  REACiOR  LFNGTH  IN  FEET* 

RK  -  1  1  COND  T IV I  TV  OF  SOLID  SP  • 

SC  —  SPFCIhIC  HFA 1  OF  SOLID  SPHERES* 

SO  -  DENSITY  OF  SOLID  SPHERES. 

SU  -  VELOCITY  OF  THt  MOVING  PFD, 

SR  -  RADIUS  OP  SOLID  SPHFRF. 

HP  —  H E  a  T  TPANSF FR  COEFFICI  E  N  i  B  F  7  w  F  F  N  GAS  A N n  S 0 LID  r-  p H F R f  3  • 

’  -  ■-  A  OF  5( iLIO  PH  °  UN  I T  1  PA  r. 

GMO  -  INITIAL  MOLAR  FLOW  RaTF  OF  GA3 * 

CPA*  <  h f -  »  C P (  —  MO L AR  oPFCIFI  »  L 

P  l  Y  OROG  FN  PESPFCTiVFL  Y  . 

VO  -  VOID  FRACTION  OF  MOVING  FED* 

TB  -  Su.ME  REFERENCE  TEMPERATURE* 

0 H 1  -  H E A  T  OF  REACTION  AT  R E F E R F N C E  T E M P E R A T U P E • 

OR  -  UN  I  VERS  I  AL  GAS  CONSTANT  (  pS  I  A*CU.  FT  *  /  LB-  MOLE  *DFG.  R  •  )  • 
h  R  Ac  i  —  P  lv  ^  ■  i  I  O  N  A  L  CON  V  t  R  S  j  ( )N  Of  R  t  A  <  T  AN  1  IN  OU  TP  U I  R  * 

TG  -  GAS  TEMPERATURE  VARIABLE. 

X  -  FRA<  i  IONAL  CONVERSION  V  A  ° I A  RLE. 
p  -  TOTAL.  pRPSSURp  VARIABLE  IN  P  FACTOR  ( *— '  J  A  )  • 

L  -  LENGTH  OP  REACTOR  FPQM  SOI  ID  INI  FT. 

ST  -  SURPASS  TFMPFRATURP  Oh  SD(  /l)  S  p  H  P  R  P  S  AT  i  . 


D  I  M  P  NS  I  ON  F  (  i  0  )  *  P  (  1  0  )  *  Q  (  1  0  )  *  R  (  10)  *  S  (  1  0  )  *  T  (10)  *  >(10) 
1  00  READ  (5*7)  i  >Z 
7  FORMAT  (  F  5 • 2 ) 

IP  (dZ.LE.u.u)  GO  To  101 
ReAd  (5*1)  SK*  o  (  *  p.  u  *  SU,  SR 
1  FORMAT  ( r  6 • 2  *  P  8  •  4  *  P8.2*  F  8 . 2  ♦  P  8  *  4  ) 

R  P  AD  (5*2)  HP*  SA*  GMO*  CPA*  Cpf*  CPC 
FORMA  I  (  !-  6*  7  *  P  8  •  2  *  3  F  *2*  F  6  *  2  ) 

R  P  A I  >  (  5  *  *  )  VO*  ■*  i )  H  *  G 1  -  *  u  R  A  r 
3  POP MAT  (  h  5  •  2  *  FT*1*  F  •  *  F  ■  *  *  F  •  F  ) 

W  R  f  T  u  (  A  *  9  )  DZ 


P  ORMAT 

(  1  H—  * 

5  X  * 

2  OH  I NI 

r i al  incrfmfnt  =  *  re.'  > 

WRITE 

(6*14) 

1  4 

FORMAT 
12HST  ) 

(  1  H-* 

8  X  * 

2hT  1  * 

9  a  *  2  fi  i  2  *  9X  *  2  HI  3  ♦  9X*  2H  i  4  ♦  ;X* 

wR  f  TP 

(  6  *  1  .8  ) 

1  5 

format 

(  I  H  * 

]  4A 

*  2  H  T  G  * 

16X*  1 HX  *  1 6 X  *  1  HP  *  15X*  1HL  ) 

<■  CALCULATION  OF  VARIOUS  CONSTANTS. 


H=  ■  /  • 


8  X  * 


•  •  » 


l\ . 


•  (  •  (  ' 


0 

0 

•  ’  • 

1  1  ^ 

0  M  '  * 

. 

— 

.  1  t 

• 

1  )  M  f  1  1  )  ^ 

(  '  • 

l  S  •  *  -t  i 

IU.  •  H  1  . 

• 

• 

.  .  *  (  i  . 

•  • 

•  • 

•  •  ^  *  S  •  ) 

.  ^ 

• 

.  .  W  1 

<  V  .  *  J 

•  • 

0  0  •  0 

0 

• 

0  •  (  '  0 

•  • 

•  • 

•0  0^0 

'ft  (  •  ^ 

•  ^ 


•  «,  •  -  I  ) 

(  ■  f  0  A  ) 

•  •  G  ,  —  M  )  I  ' 

(  “  )  »  ) 

•  «  .  •  1  !  I 

0  \ 


AH=Sk  /  l  SC*Sl)*SU*H*r-i ) 
M  7 =  p  J#  ( )*m*HF / SK 
GAat99?.0/(  1P4.0+RT  ) 
f  A  =  9.  )*GA-42. 
f R  =  -9  • 0*GA+46  • 

PC  =  2 2 P • 0*6A-1 410.0 
F  D  =  0  A  *  R  T 

c  a  =  hf*.sa/gmd 
r pM=-r pa  +  2 . Q*r pc +  r K - 

(P=i .0/ (GMO*GR ) 
XA=SOPT (2.0) 

U  A = ( X A- 1 .G)/2.  ) 
UB=(2.0-XA)/2.0 
U  C  =  —  X  A  /  2 . 0 
IJ I )  =  (  2  •  0  +  XA  )  /2  •  C 
n  L  =  0 . 0  0 


[  NHUT  OF  I  N  T  T  1  A| 

da  r A  . 

P  F  A  0  (3,4)  (  V (  I  ) 

* 

- i 

II 

s 

JO 

4 

FORMAT  ( F  9 • 3 ) 

.  s  . 

1  1 

IF  (C0UNT.GT.iu. 

0 )  GO  TO 

1  3 

Z  = DZ 

GO  TO  18 

1  3 

IF  ( COUNT  *  GT • 2  . 
Z=OZ+07 

GO  TO  18 

0)  GO  TO 

17 

1  7 

Z  =  3 . 0*nz 

1  8 

DO  IP  1=1,8 

1  Q 

T  (  T  )  =  V  (  T  ) 

P  C  =  1  .  0 

C  .SOLUTION  OF  PRORLFM  BY  RUNGF-KUT  T  A-G I LL  METHOD. 

2  0  F (  1 ) =  A  P* (  —  26.  )*  i  (1  ) +2 “ • 0*  T  (  2  )  -  (  ^  )  )  /  R  ,  0 
F  (  2  )  =  A  p  *  (  -  '  11  )  —  I  2  •  0  *  T  (  2  )  + 1  3  •  0  *  T  (  3  )  )  /  9  •  0 
r  ( 3  )  = AP* ( 2 • 0*T ( 1 ) +3 •  *T (2 ) -24 • 0*T ( 3 ) +1 9.  T ( 4 )  )  / 1 5 . 0 
F ( 4 ) s AP* ( 2 • 0*T { 2 ) +9 . 0*  ) - 36 • 0* T ( 4 ) +2  . 0*T (  )  )  / 2 i • 

>  =  (  9. G*T ( 3 ) -  0.0*  (4 ) -  •  >  )  +67*7 ( 6 )  ) / ( i 64 • 0  +  T ) 

T  K  = ( 12.43-22100.0/7 (6)  )/0.434294 
T  R  =  3  6  u  0 . 0  *  F  X  P  (  T  <  ) 

F ( 5 ) =  A  P* ( 99  2 . 0*  T  S— 1 4 1  •  •  ■  •  -  (  ) +46  •  (  4  )  — 4  2  •  (  )  ) / 2 

F(  7  )  =  C  R  *  T  (8)*(  1  .0— T  l  7  )  )*TP/(  T  (M*(  1  •  >  +  2  •  1  *  T  (  )  )  ) 


F  A 

n 

(  (  T 

(  6  ) 

-T 

R 

)  mm  +  i  ) H T  )  *F  ( 

7  ) 

F  ( 

6 

)  =  ( 

C  A* 

(  T 

—  T 

(  6  )  ) 

—  h 

A  )  /  (  C 

P  A  +  T  (  /  )  *  C  P  3 

F  ( 

8 

)  =- 

•  i 

1* 

( 

I  . 

0+ 2  . 

U* 

T  (  7  )  ) 

*  T  (  6  )  /  I  (  8  ) 

I  F 

(  PC 

. 

•  i 

• 

0  ) 

GO 

T  0 

2  3 

DU 

22 

1=1 

,  d 

P  ( 

I 

)  =Z*F  ( 

1  ) 

22 

I  ( 

I 

)  =  Y 

l  1  ) 

+  p 

( 

1  ) 

/  2  •  U 

PC 

= 

PC  + 

. 

GO 

TO 

20 

2  3 

1  F 

(  PC 

.  GT 

.  2 

• 

u  ) 

Go 

To 

26 

no 

28 

f  =  I. 

,  F 

\  (  1  ► 


. 

*  .  - 

L'  . 

.  *  r  -  . 

(  OrHO  '^)  \  .  r  = 

(  .  )  or 

.  s  \  (  •  f  -  x  )  = 

.  \  \  (  x  - 

.  \  X--01 

AM  y  +•  .<•)=<" 

i.O*  |#1 

l  .  r  r  .  (  )  V  )  (  .  ) 

(  •  #  i  ) 

(  •  •  • 

f  .  m  .  r  .  1/ 

Vi*  .*■  = 

-  •  ’  -  i  f 

(  t  >  y  *  (  »  )  > 

#-\{(  -(•>.)  .  •*■  <  )  *  ‘  a  -A  <,  —  )  A  *  ( 

(  (  *•  )  7*i  •  ►  f  +  <  v  J 

.  4-(u)  i  •  -  ( 

.  -  )  (  (  ^  )  (  )'■"■•  *■  t  •  •  +■  •  '  ~ 

•ji  +1  #  \  (  (  *  I  \  •  v  —  ‘  1  • 

I  >  •  )MX  1*  .  = 

i  .  v  \  -  <  u  )  .  4  i\  +  <  i  •*  .  *  1  -  .  ‘ 

.  4-  .  I  )  *-  (  )  I  )  \  +■(('  I  ~ 

* ■  i.  »  M  *  (  —  (  A  I  )  )  •=  > 

I  <  •  )•  .  V  +  .  )  >'  1  v  • 


.MU  iHf  I  |)»r(  |  J  1 

.  «  - 


P  c 

= 

PC 

4- 

1  . 

(j  ( 

I 

)  = 

2  *h 

76 

1  ( 

I 

)  = 

Y 

(  1 

GO 

1  0 

70 

26 

1  F 

( p( 

.G 

no 

7* 

1  = 

K  l 

1 

)  = 

Z  *r 

2  8 

T  ( 

1 

)  = 

r 

(  ) 

PC 

= 

PC 

-i- 

1  . 

GO 

(0 

2  0 

8  u 

DO 

31 

I  = 

S  l 

I 

)  = 

Z*t- 

31 

Y  ( 

I 

)  = 

V 

(  ( 

c 


)  +  UA*P  (  1  )  +  lJ H *Q  (  J  ) 


3  .  ) 


)  +U|)*R  t  i  ) 


3  3 
8 

10] 

$  h  N  T  P  Y 

0 .  u  2 

1  •  8  0 
#7.6 


)  +  (  p  (  I  )  +5  (  I  )  )  /6  •  0+  (  Uh*<j  (  i  )  4- u i  >*p  (  I  >  )  /  3  •  0 

i ) L  =i  >L.  +  Z 
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1  =  1 

7  2  W  P I T  P  ( A  *  4 )  I  j  M F 

4  FORMA  1  (  1  H  U  <  4HAI  ^  F  H  .  4  ,  AH  MIN.) 

W  R  J  f  F  (  A  y  7  ) 

•  1  ( *  H  ’  RX»2HTG>] 2X  *  2HT 1  *  7X  *  ^ H 7 2  *  7  X  >  2H 1  -  .  X , 2 H  4  *  X  *  2 1 
calculation  op  constants  FOR  ThF  matrix. 

a u  A N  =  n|  * 3  .  u*  (  l  •  0  —  V  )  *  F  I  * H P  /  (  R*V  ^\/F*rp*i')p  ) 

» N  =  Ft  *s< /  (  v  f  *  r  s  *  n  s  *  o  R  *  n  R  ) 

CN=a j.  *hd*hf /S< 

OM  =  QO7#0/ (  1 p  4 . 0  +  r n > 

DB=nT*PM/2.0 
F (  1  )  =OR/P . 0 
F (  ? )=  —  DR/7.A0 
F (  3 ) =  —  I )  R / 10.50 
F ( 4 ) =42 .0-9 . 0*DN 
F  (  5  )  =  5  0  •  u 
E ( A ) =-9.0 
4  )  A  (  ]  )  =n  .  0 

A  (  ?  )  =!)R/9 . 0 

A  (  3  )  =  —  D  P  /  A  .  0 

A  (  4 ) =-Q . 0*OP / ? 1  . 0 

A  (  A  )  = AO  . 0  *  RM  —  4  A  0 . 0 

A  (  A  ) =  —  225  .  J 
B(  1  1  =  1.0+13. 0*1  )R  /  4  .  O 
R  (  2 ) = 1  . 0  +  4 . 0*OR /  3  .  J 
R  (  3 ) =1 . 0+R.0*OR/ A.  0 
R ( 4 ) =1  . 0+1 2 . 0*OH / 7 . 0 
R  (  5 )=270.0/DR+]41 0.0-225.  '*DN 
R  (  6  )  =  T  P4.  J  +  2 . 0*  (  1  8  4.  o  +  (  j\|  )  /  AN 
C  (  1  )  =  —  27.  ()*OB  /  8  •  0 
c(?)=-n.  o*r>B/Q.  o 
r  (  3  )  =  - 1  Q .  0  *  O  R  /  1  a  .  o 
r ( 4 ) =-? a . 0 * O R / 7  l  . 0 
r (  a  )  =-nN*rN 
4  1  C  (  A  )  =  0  •  0  0 

n (  1  ) =  T (  l  *  T  ) * ( 2 . 0-M (  1  )  )  +  T  t  2  *  (  )  *  (  ~r  <  1  )  ) +  T ( 3 . i  ) * ( - F (  1  )) 

D( 2)=-A(2)*T(1 *l)+(2. 0-3(2) )*T (2*1)-  ( 2 ) *T ( 3  *  I  ) 

D(  3  )=-P  (2  )*  I  (  1  *  I  )  —  A  (  4  )  *  i  (2*1  )  +  (  2  •  0  — B  (  3  )  )  *  i  (3*1)—'.  (  -•  )  *  i  (  u  *  I  ) 

I)  (  4  )  =-P  (  3  )  *T  (2*i  )  -a  (  4  )  *  T  (3*1  )  +  (  2  •  u-B  (4))-*T(4*I)-f(4)*T(r  *  I  ) 
f)  (  5  i  (  4  )  *T  (  3  *  1  )  —  A  )*|  (4*1  )  +  (270.G/  +22  •  4  .  0  )  *  (  5  *  ) 

I)  (  5  )  =n  (  A  )  +TG  (  I  )  * l ) N * C N 

n ( A ) =9. 0* I  ( 3  *  T  ) -50.  )*T(4*T)+22A.n*l  ( u .  I  ) 
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P  C 1 =P (  1  )/(  (2) 

R (  1  )  =R (  1  ) -A ( 2 ) *P  C  J 
(  (  1  )=C(1  ) —  R ( 2 ) * P (  1 
n  (  l  )  =  n  (  1  )  -i)  (  2  )  *F  d 
FT2=F(A)/F(4) 
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1  5  0001  REGtNFRAT  IV  E.  t.  i )  nFAi  EXCHANGER  —  Fn  COUPLING* 
004 

FIXED  GO 

EXCHAN  NOL  I  ST  *NOf)EC< 


D  t  5  I  G  N  0  E  R  u  G  E  N  r  R  A  T  I  V  F 
IN  SOLID  SPHERES  - 


n •  >  H-  A  i  f  a  HA N G ?■  w  w  I  i  H  L  u  »<  i  H r  p m A L 
RY  FHE  MFTHOD  OF  dECOUPL  I  NG# 


d  I 


r—  I  • 


ALL  I NPUl  DATA  \R  F  IN  LH  *  • j  -  »  F I *  B TU  *  DFG  E  *  * 


FOR  C  0  N  v  F  R  G  F  im  C  Y  uE  I  H  E  NUMERICAL  INTEGRATIONS  7  Mr 
DIVIDED  jNTO  XJ  equal  sections, 

INTEGRALS  ARE  EVALUATED  BY  SIMPSON  METHOD  WITH  ERROR 
D  I  Ft- 1  REN  i  iAL  FQG'Ai  IONS  ARE  INTEGRATED  RY  RONGE  —  Ku  T  T  A 


-  S  SpF  <  • 

RAi'IUS  HAS  r  ■ 

CORRECT  I  ON • 

-  H  0 1  • 


Gr  -  HASS  FLOW  RATE  OF  GAS, 

EL  -  LENGTH  OF  Thf  REGENERATIVE  RED, 

To  -  INITIAL  TEMP FRA TURF  OF  SOLID  SPHERES, 

V  -  FRACTIONAL  VOTD  OF  THE  RFGENpRATIVF  RED. 

A R  -  RADjUS  OF  Thf  SOLID  S P H F P F , 

DS  -  DENSITY  OF  SOLID  SPHFPFS, 

CS  -  SPg<  I  FTC  HEAT  op'  Thf  SOLID  SPH  RES# 

hi-  -  f*f  A T  TRANSFFR  COEFFICIENT  BETWEEN  GAS  AND  SOLID  SPHERES, 

S K  -  THERMAL  CONDUCT  I VI  TV  OF  SOLID  SPHERES. 
cp  -  Specific  heat  of  gas  at  mean  tfmhfraturf, 

DENSITY  OF  GAS  AT  MEAN  FEMPFRATl  . 

EC  r  —  part  of  a  cycle  required  IN  THE  program  in  minutes# 

CT  -  TriF  TIME  E'ER  CYCLE  IN  MINUTES. 

Tg(  1  )  -  INL^T  GAS  TFM  P  F  R  A  T  UR  E • 

ES  -  FRACTIONAL  FRRGR  ALLOWFd  IN  THF  ENERGY  BALANCE  BETWEEN  GAo 
AND  SOLID  SPHERES, 

FNl  -  NO,  OF  i  MCRFMFN TS  PrR  t  YFL*-, 

NX  -  NO.  OF  daTa  iNPul  F 0 R  ERF(X). 

N  3  _  21  GIVEN  d i V i SIGNS  A  L  0 N G  THE  R A  f )  I U S  0 E  S P  H F R E • 

im 4  -  NO.  OF*  i  N C R  E  E N  T S  A  L 0 N G  THE  ;  ^  A  '  7  0  .’•>  . 
hRF i a)  -  ERROR  F  U  N C  T I  0 N  AT  A , 

Tu  -  GAS  TEMPERATURE# 

T(J)  -  SOLI',)  TEMPERATURE  AT  THE  JTH  SECTION* 

TS  -  SURFACE  TEMPERATURE  OF  SOLID  SphERi • 


OUTpUl  DATA  ARE  IN  7  COLUMNS  OF  GAS  TEMPERATURE  AND  SOLID 
Ur  ID  T7*  •  :  *  r  1 F ,  7  1  9  AND 


DIMENSION  T G ( ? 0 )  *  T ( x o  *  X 1  )  *  TT ( X  M  »  P ( X  I  >  *  F x A ( 2 1  )  *  F X (  ( 2 1  ) 

D  I  MENS  j  ON  BEk(21)*  PA( 21*21)*  p6(21*21)*  a  125)*  EREi2“)*  »'  (  2  1  ) 
1)1  Mr  NS  I  ON  S(a)*  PE  (21  *21) 

INPUT  DATA. 

READ  (5*1)  bF  *  El*  Tu*  V*  AR *  Do*  Co*  HF 
1  FUKMAi  (  P  8 . 1  *  r  6  ,  x  *  F  7 • 1  *  F  7 , 5  *  F T ,  :  *  F  7 . 1  *  F 7 • 3  *  F  .  2  ) 

R  E  A  D  (5*X)  S K ♦  Ch  *  D  F  *  Cl*  F C 7 
X  FukMAi  (  F  8 • 3  *  F  I • 3  *  3  F  /  •  2  ) 

READ  ( s *  i  )  T  b (  1  >  >  F  S  *  F  N 1  *  N  2  *  N  7  *  N  a 
3  F  UKMA i  ( F  , 1  *  F  8  •  4  *  F 8  «  2  *  7  I  a  ) 

W  R  f  T  F  (  A  *  H  4  ) 

Ha  FORMAT  (  1  H- *  luX*  X8HDATA  GIVFN  EELo'W  M~r  E  -  ) 

W  R  T  I  F  (  A  *  «  5  ) 


» 

• 

• 

• 

# 

J 

.  • 


9  f  J  • 


•  < 


•  ' 


.  •  J  ~  •(Is#  i  •  ( 

f  (  s  ) 


*  ! 

M  •  ( 

• 

9 

•  < 

•  <  « 

(  i  • 

• 

•  • 

•  • 

(  • 

(  • 

•  • 

• 

• 

#  •  • 

• 

• 

•  • 

•  . 

.  r 

• 

»  4 

I  •  A  , 


h  O  R  H  A  I  (  1  H  0  »  a  9  z  H  i  G  9  9  a  9  <+  H  I  (3)9  6  A  *  4  H 1  (  7  )  9  o  A  9 

1  8  PI  I  (  1  6  )  9  b  A  9  b  H  1  (  I  9  )  9  A  a  9  2  H  1  b  ) 

A I  ■  L  '  i  ON  Oh-  VAR  i  OUS  C  ON  i  A  .  i  . 

FN3=N3-1 
DR  = AR/FN3 
F  N4  =  N4 
1)1=1  .  J/FN4 
I  1  vih  =  U  .  Uo 
T  K  =  S K /  (  !.)S*Cb  ) 

VF =  GF/ ( DF*V ) 

P  I  =  3  •  1  4  1  d  9  3 
dT  =  C  T / ( 60 • 0*F N 1  ) 

R  t  A  0  (  6  ♦  4  )  i  a  (  j  )  9  F  k  F  (  J  i  *  J  =  1  9  N  2  ) 

F  0  R  v!  A  7  (t-7.^9  ^  A  9  F  8  •  6  ) 

A  1  =  3.  J*oZ*F  L*SF  ■*  i  I  ♦  0  —  v  >  /  A  < 

A  6  =  SQR  7  iDT/(Pl*TK)  ) 

A 3=-oT/AR+Z*0*A6 

A4  =  -dT*I)T  /  1  ;.u*ak  )  +4  •  0*0  7-*  Ab  /  •  0 

A  5  =  bOR  1  l D  f*P 1 *T K  ) 

A  (  =  ]  • J /  l 4 • 0* T’K*D  i  ) 

A  8  =  1  •  u  /  (  ^  •  U*S(JK  T  (  I  <  *  D 1  )  ) 

1  =  d  I  *  V  F  *  V  *  C  P*  D  F 

.  .  -  F* t L /  (  AR*V-»CR*uf 

H 3  =  1  •  0 /  (i£:»G*A')*7  <  *  D  7  ) 

H4  =  iZ»  U*DK*hF/bK. 

DO  9  i  =  I  9  IM 4 

D()  Q  J  =  I  9  N  8 

T  (  T  9  J  )  =  !  O 

DO  10  1=1  9 N3 

F 1=1-1 

R ( | ) =F I *DR 

DO  11  1 =  j  9 N 8 

t  A  A  (  i  )  =  0  •  0  0 

F  A  C  (  i  )  =  0  •  U  0 

DO  11  J  =  J.  9  IN  3 

DA  l  I  9  J  )  =0 . UU 

PF  l  1  9  J  )  — u*0u 

Phi  1  9  J ) =  j • 0  u 

DC)  18  I  =  1  9  N  8 

UA  =  R (  1  ) *R (  l  ) *A  7 

UC= l AR-R ( i ) ) **?*A7 

fh  ( UA • O  f  . SO . 0  )  GO  TO  17 

h  A  A (  1  )  =  ]  •u/h  a  p  i  u  A  ) 

it-  i  Oc  .GT.bO.G)  bO)  TO  1  ' 

F  A  c.  I  1  )  =  J  •  0  /  t-  A  P  l  U  (  ) 

CON  7  1  IM u  P 

DO  2u  1  =  i 9  N3 

DO  2  0  J  =  1 9  N 3 

U  A  = ( R (  1  ) - R ( J  )  ) ** 2  *  A  7 

U  P  =  (  R  l  I  )  +  R  (  J  )  )  *  *  2  *  A  1 

U F  =  (  - 2  •  o*  Ak  +  R  (  1  )  +R  (  j  )  )**</* a  7 

1  1-  (UA.GI.  30.0)  00  TO  2 

PA  (  1  9  J  )  =1  .  O/'F  Xp  (  UA  ) 


1 

5  H T (  1  1  )  9  7A9 


* 

.  =  1  * 

•  *  • 

.  )*•£+!-  .  - 

•  4*  J  s 


• 

• 

• 

• 

•  r 

♦ 

• 

• 

•  ' 

# 

• 

\  l  • 

(L» 

• 

• 

• 

• 

• 

(  '  * 

i  —  I  U  4 


2  2 

I  F 

(  U  r 

•  GT 

•  50* 

)  ) 

GO 

T  n 

1  0 

2  4 

h6  ( 

i  *  J 

)  =  1 

•  0  / 1 X  P  ( 

) 

2  a 

I  F 

(  Ur 

•  (t  ! 

•  9  C  • 

1  ) 

GO 

TO 

20 

HF  ( 

I  *  J 

>  =  i 

.  u  / 1-  XP  (  UF 

) 

20 

CONTI NUF 

NO 

1  16 

1  = 

1  *  N  3 

X  j  = 

( AR 

-R  ( 

I  )  )  *  A  8 

119 

DO 

120 

J  = 

1  *  N  2 

I  F 

(  X  J 

.  LT 

•  X  (  J  ) 

) 

GO 

TO 

121 

I  F 

(  Xj 

.  EQ 

•  X  (  J  ) 

) 

GO 

TO 

12  2 

12  0 

CON 

T  INUF 

XER 

=  C.  . 

0  0 

00 

TO 

117 

122 

XER 

=  1  . 

0-E 

R  F  (  J  ) 

GO 

TO 

1  1  7 

121 

NX  = 

1  X  j 

-X ( J- 1 ) ) 

/  6 

•  0 

N 0  12  3  JA  =  1  >  7 
FJA=JA— 1 

X  X  =  (  X  (  J  —  t  )  +F  J A-*-nX  )  **2 
12  3  V  (  J  A  )  =  1  •  vj  /  F  X  P  (  X  X  ) 

XE R  =  41 • 0* ( V ( 1 ) +Y ( 7 ) ) +5 1 6 • 3* ( Y ( 2 ) +  i 6 ) ) •  .  (  ( 3  ♦  (  5  )  ) 

XER=(XFR+2  72  •  0*Y  (  4  )  )  *  i )  X  2  •  )  /  (  l40*0*SQPT  (PI  )  ) 

X  E  K  =  •  0  —  r  R  e  (  J  —  i  )  —  X  E  R 
117  BER(I)=XER 
116  CONTINUE 

DO  7  3  1=1 » N  3 

FXA (  I  )  =R (  I  ) *R (  I  ) *EXA (  I  ) 

1)0  7  3  J  =  1  ,N3 

73  PA(  I  *  J ) = ( PA (  I  *  J ) -PR (  I • J ) -PF (  I  ♦  J  )  ) *R ( J ) 

102  I A= 1 


8  1 
C 

1  ol 


7  I  M  F  =  T  I  M r"  +  0 7-^6  0  •  0 
wK  I  it-  (6*81)  1  i  M  t 

F  0  R  M  A  1  (  1  H  0  *  3X9  3FiAl  *  F  9  •  6  *  6H  MIN*  ) 

ITERATION  FOR  ENTHALPY  BALAN<  r  E  TWEFN  THE  T w’O  PHASES. 

P  =  T  (  I  A  ♦  N  3  ) 

f  S * B4*TG  (j  I A  )  + 1  8  •  0*T  (  I  A  *  N  3- 1  ) -9  •  0*  (  1  *  N 3-2  )  ■*■  7.  •  •>  (  1  *  -  j 
T  S  =  T  S / (  1  *0  +  64) 


0= ( TS-P ) / NT 
DO  A  o  T  =  I  1  *  N  3 

30  Y (  T  ) =  R (  r)*T(IA*I)*(-RFP(I)  ) 

A  2  1  =  NR  *  2 . 0*  (  Y  (  2  1  )  +  2  •  o*  (  Y  (1  3  )  +  Y  (  ]  7  )  )  +  4  . 
A22=2*0*(  Y (II  ) +  Y ( 1 3 )  +  i ( i  9  )  +  Y (  "  7 ) + V (  '  0 ) 
A  2  2  =  F>R*  ( A  2  2+4*  o*  (  i  (  1  8  )  +  Y  (  2  0  )  )  +  Y  (  2  J  )  )  /  3 
A23=(A22-M2J  )  / 1  3  •  0 
A 2 = ( A22+A2 3 ) / i AR*TK ) 

l^u  Ui=Al*'(A2+A3*F  +  A4*'0) 

\  A  =  T  G  (  1  A  ) 


DO 

40 

19  =  1 

*  4 

S  (  K 

)  = 

B2*NZ* ( P+Q*NT / 2 

1  F 

(  K 

•  • 

2  ) 

GO  IU  42 

I  A  = 

1G 

(  F  A  )  +  S  (  K  )  *  0  .  ^  u 

GO 

TO 

4  0 

[  F 

(  K 

•  (1  I  • 

3  ) 

GO  IO  4  *■> 

0*  (  y'  (  1  ]  )  +  V  (  ]  k  )  +  V  (  1  q  )  )  )  /  3 
) +4 • 0* (  Y  (  1 2  ) +  Y  (  i  4  ) +  V  (  In) 

•  o 


.  u 

) 


•  ) 

. 


1 1  . 

« 

• 

.  L- 

. 


. 

-  •  - 

t 

(  A  X  i  H  A  X  ^  • 

•  1  £  = 

-  —  • 

*  r  s  } 

- M *  fsL  r' 

f  i  •  i  }  --\t  *M,k  « < (  •T)A,m)  =  (C*T) 


•  • 


• 

• 

9 

t 

• 

• 

• 

•  * 

*  • 

A# 

•  * 

.  ) 

1  -  1  0  5 


I  A  =  I  G  (  T  A  )  +  S  (  X  ) 

40  rONTlNUF 

45  T G ( T A+ 1  )  =  TG(IA)  +  (S(  1  )+S(4)+2.0*(i  ( ,  ) +  S ( 3 )  )  ) / 6 • 0 

Q  ?  =  (  T  G  (  I  A  )  - 1  G  (  1  A  + 1  )  )  *B  ] 

OS  T  =0*0  T  4-  p 
T  h  S  =  A  H  S  (  O  1 -02  ) 

T(r  S=  f  F S *2  •  0 /  (  0 !  +Q?  ) 
if-  (TFS.LF.FS)  GO  TO  AQ 
54  Q= (  (01 +02*2 . 0 ) / ( 3 • 0 *  A ]  ) - A  2 - A  9  *  P ) /A4 

GO  TO  iOu 
3  u  00  1=2  *  N3 

00  5  A  J=1  *  iM 3 

A A  Y(J)  =  r(TA*J) *  P  A (  |  . J ) 

W1  =4.  )*Y  (  3  ) 

HO  A  P  K  =  5 .  1  7 ,4 

a  p  W] =W1+2.0*Y (K )+4.0*Y (Y  +  2 ) 

A  2 1  =0P*2 . U* (  Y  (  7  )  +  Y  (  2  1  ) +  A  I  )  / 3 • 0 
W2=4.0*Y  (  2  ) 

00  3  0  K  =  :>  *  I  9  *  3 

39  W  2  =  W2  +  2  •  0*  Y  (  <  )  +4. 0*  Y  (  K  +  1  ) 

A22=0R*l Y ( 1  ) +  Y (21  ) +  a  2 ) / 3 • 0 
A2  3=  (A22-A21  )  /  1  3 • 0 
T  A  =  (  A  2  2  +  A  231/(2.  U*R  l  i  )^A3  ) 

T  A  1 =AR*P*RFR (  I  ) / K (  I  ) 

TR=(I)T+(AR-R(  1  )  )  **  2  A(  T  K  +  T  <  )  )  *  R  F  R  (  ]  ) 

I  r  =  -  (  A  P  -R  (  I  )  )  *  A  A  *  F  X  r  (  T) 

7  r  i  =  i  i r  +  r  r ) *  a  r  *0 / p  (  i  ) 

a  A  TT  (  T  )  =  T  A+  I  A1  +  TG7 


P  3 
C 


T  l  (1  )  =  2  P  0 . 0 

w  R  I  I  F  (  A  *  H  3  )  T  G (  i A  + i  )  *  T  T (  3 )  .  ,1(7)*  T  i  (  li  )  * 
FORMAT  (  7  H  *  IX*  F7U.3*  IX*  AF70.2) 

FND  OF  ONF  RFAGTOR  I NCREMHMT  CALCULATION. 


I  F 

(  I  A .GF • N4 )  GO 

TO  61 

o 

o 

A  4  1  =  1  *  N  3 

A 4  T  (  I  A  *  T  )  =7  1  (  f  ) 

7  A  = 

T  A  + 1 

GO 

TO  1  07 

AT  T  F 

(  7  TMF.LT.FGT  ) 

GO  In  702 

WP  T 

IF  (  A  *  P  !  ) 

p  7  FORMA  1  (  1  H— ,  TAX* 

2  1  H  F  N  D 

OF 

f  n  n 

• 

sfm  try 

F  XG HAN 

2000 . U 

4  .  A  P  2  P 1  >  .  0 

.34  . 

0  3  0 

O.lQO 

0.2^0  23.00 

a  .  no 

35  ).0 

0.0010  4  p . 0  0 

2  1  2  J 

7  3 

u  • : )  o  u 

• 

U  .  1  u  u 

0.11 24 A  3 

.  ? 

0.222703 

.  3  )( 

0  •  3  2  P  A  2  7 

1.400 

0 .428302 

n  .  c;  o  n 

0  .  A  2  0  A  0  0 

0  .  A  U  ( ) 

0  .  A  0  -2  p  c;  A 

1,7.10 

0  .  A  7  7  p  0  1 

'ALCULATION.  ) 


(7Q4.  A  0.191 


l  i  (  i  3  )  * 


2 ' '  .  0  0 


1 7  ( 1 9  )  *  ns  t 


*  (  * 

• 1 


. 

- ^ ' s=  t 

(l*|)  q*  {  I.  •  AT  )  I 

•  '  ,  f  # 

•  A ♦  (  >  >  V  *  .  4 

(  S  J  V-v  .  4i  - 

•  f 

. 

. 


' 

(Ml 

•  ( 

. 

*  v  • 

f  • 

• 

* 

*  ^  ) 

*  — 

o  '  . 

1  -1  06 


0  .  ft  00 

0.7421 01 

0 .  o  on 

0.7QAO0P 

i  .  00  n 

0 . p  427  01 

i  .  7  0  0 

O.oi 0214 

1 . 400 

0 . 9  5  ?  2  R  5 

1  .  A  0  u 

0.976348 

1  .80  0 

n.ORQOQ] 

2  .  (.'  JO 

0.995322 

2.200 

0 .  •  8  1  3  7 

2.4.) 

0.9993  1  1 

2  .  A  U  0 

0.99 9 744 

2  .  ft  0 

il  OQano  5 

2.0  0  0 

0.QQOQ7P 

(SAMPLE  OUTPUT) 


150001  RFGFNFR AT  I VF  RFO  Hh  a  l  FXTHANGFR  -  DECOUPLING. 


AT  Q.POOPPO  M  I  N|  . 


3  0  p 

• 

063 

28  1  .67 

783.44 

7  0  8  .  r  1 

2  °  8 . 8  6 

616.17 

324.82 

8y  7 

• 

5^2 

781.11 

287.41 

2  8  6.  ] 4 

294.1 4 

307.60 

-  1  . 

3  1  8 

• 

3  96 

280.76 

281 ,68 

284.42 

290.64 

301.35 

3  0  6  .  6 

31  o 

• 

707 

280.51 

281.17 

2  «  3 . 1  7 

287 .82 

296.46 

302.61 

3  u  4 

• 

326 

280. 34 

280.81 

282.26 

285,77 

292.60 

207.69 

299 

• 

107 

280.23 

280.56 

281.61 

2  8  4  •  2  4 

269.68 

2  9  *64 

294 

• 

892 

280.16 

280.38 

2  8  1 .14 

2  8  3.  1.0 

. 

290 . 62 

29] 

• 

3  28 

280 . ] 0 

280. 26 

2  80.8  1 

282 .26 

2  8  6 . 46 

288.05 

288 

• 

8  6  1 

2  8  0  .  o  6 

280.18 

2  8  0.6  7 

281 , hn 

284.06 

286.14 

286 

• 

777 

280.04 

280.12 

2  8  4 . 4  o 

281.19 

268. 04 

,  84.6! 

28  6 

• 

1  47 

280.02 

2  8  0 . 0  8 

260.28 

280.86 

2«  6. 6 1 

28  3 

• 

8Q7 

780.01 

2  «  r'  •  0  6 

7  p  0  .  1  9 

260.62 

■  i  .66 

282.63 

287 

• 

467 

270 . QQ 

2  7  0.09 

2 7 o.qq 

2  8  0 . 0  0 

280.34 

28  1  .61 

A  1 

083333  M [ N . 

3  4  3  .  /0  7 

3  0 1  •  f 

3 0 5.27 

311.78 

320. 6 1 

3  6 ]  .24 

• 

3  •  7.2 2  - ) 

2  9  7.32 

300.2 7 

306.08 

314.41 

324.53 

32  -  •  ? 8 

330.791 

293.63 

296. 1 7 

8  0  1  .26 

-  1  -  . 

3)6.8  3 

:•  2  8  ,hf- 

324.  i 

290.68 

2  -  . 

2  97 • 2  3 

30  3.  ■ 

312.72 

3  j  7. 82 

288.33 

2  O  )  .  i  - 

293.89 

299. 75 

307.7] 

812.47 

3  1  3.464 

286.46 

2  -  7 . 

?  9  1  .  - 

296 . 2O 

303. 80 

3  0  >  . 6  - 

R  .  6  )  7 

2  8  6 . 0  U 

2  6  6.2-8 

266.64 

208.21 

2  '8  0.46 

30  8.  3  6 

304. 2 7 o 

268.66 

2  «4 . 6  6 

. 

2o0. 72 

■  • 

79.6  3 

3  U  0 . 4  6  2 

267 . 06 

28  3  . 

2  6  c  .  6  6 

•  6  4 

3.31 

2  6  6.  ;  >- 

•  . 

2  8  7 . 2  6 

?  .  7  2 

2  6  4.  9 

6  •  9 

)  • 

• 

281.72 

. 

2  6  -<.44 

2  6  5  .  6 

2  6  8  . -u 

7'  •  2  1 

. 

283.31 

/  -  1  •  7 

. 

34*4  3 

• 

2  8  .24 

28  7.422 
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. 

2  7  . 

280 • 0 I 

281.07 
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1  5UU01  ACF  I  Y|  FNE  PRODUCTION  MY  i )E< HJPL  I  NG . 

QUA 

LFUNG  GO 

E  [  NAL  NUL  1ST  *  NODE  (.< 

IHlS  PROGRAM  is  de  v E  LUpe  O  IN  ACCORDING  TO  THE  R  F  GF  iN  r  7  A  T  J  V  i 

HEAT  EX<  MR  BY  DECOUPLING  TECHNIQUE*  AND  IS  SPECIALLY 

APPl  I  ED  FO  THE  PROl)U<  I  ION  OF  ACETYLENE  FROM  METHANE. 

A  |  |  T  l\!  P u  I  i  >  A  I  A  ARE  I  i\i  I  h  9  HR  *  El*  H  :  1  -  »  I ) E.  G  t-  *  1  <\i  [  i-  •  »  . 
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THFRMAL  CONDUf T { V 1 Tt  OF  SOLID  SPHERES. 

CS  -  Spec TFIC  heat  OF  Solid  Sphepes. 
ns  -  de.NStTY  of  sol  id  Spheres. 

A '■?  -  RADIUS  OF  SOL  ID  SPHERE  IN  INCHES. 

HE  -  HEAT  TRANSFER  COEEEK  IENT  BETWEEN  GAS  AND  SOLID  SPHERES. 

AS  -  TOTAL  SURE  ACE  aREA  OF  SOLID  S  p  H  t  K  E  S  PER  UN i T  VOLUME  SPACE. 
V  -  FRACTIONAL  VOID  UE  T hE  RFAClOR  BED. 

GF  -  MASS  FLOW  R  A  1 E  OF  GAS. 

r d  i  ,  <;  p  / ,  rps  -  L A R  S P E (  T  F  I  <  h  F  a  T  OP  M E  T H A N E  *  A C E  T  I  * 

AMD  HYDROGEN  RESPECTIVELY. 

FI  -  LENGTH  OF  THE  REA (  LOR  BED. 

TR  -  RASE  TEMPERATURE . 

HO  -  HEAT  OF  OF  A  FT  TON  AT  JFMPE MATURE  I  M, 

T GO  -  INI  FT  GAS  Temper a  Tup- . 
fO  -  TNJTTAI  SOI  ID  TEMPERATURE. 

PO  -  INI  FI  IOIAI  PRESSURE. 

XO  -  INLET  FRACTIONAL  CONVERSION  OF  Er Yd. 
dT  -  The  Time  per  ( yCLe  In  minute: . 

FT  i  -  pari  of  A  D  yr  I  e  Pi-ouiRFD  IN  THE  PROGRAM  IN  MINUTES. 

ES  —  FRACTIONAL  ERROR  ALLOWED  I N  I  he.  tNERGY  RaLAN(  F  BE  I  w  Fn  da 
and  SOLID  spheres. 

N  2  -  NO.  OF  DATA  input  fop  ERe(X). 

N a  -  NO.  op  DIVISIONS  IN  The  Radius  vaRIAhle. 

N  a  -  NO.  OF  TNCPFMFNTS  A I  ON  6  The  re'A<  OR  Hu. 


H  J  MENS  T  ON  1 

(  ?  F  ,  1 

7  )  ,  T  I  (  1  7  )  *  R  (  l  7  I 

% 

F  X  A  (  1  7  )  . 

FXY 

(  1 

n I  MENS  T  ON  PA  (1  7, 

17),  P P (  !  !  *  1  7 )  ,  PF (  1  7  *  1  7 )  * 

X 

I  2  f 

)  * 

DI  MFNS TON  T  Y  (  R  )  * 

,S  (  P  ,  A  )  *  E  (  M 

1 

FORMA  1 

l?Xi 

E  A  .  S 

,  Eft  .  A ,  E  A . 2  *  "A, 

2 

9  E  1  .  /  *  E 

7 

.  2  * 

F 

2 

format 

(2X9 

E  A  . 

,  E  ft  •  'J  *  E  *  •  P  *  F  e  • 

p 

•  F  (  .  2  9  e 

p 

.  I  • 

F 

P 

FORMA  f 

(  2  A  » 

E  7 . 2 

*  E  4  •  /  9  E  8  •  2  »  Eh. 

z 

9  E  7  •  /  9  F 

.  A  ) 

A 

FORMA  1 

(IX* 

E  7 .4 

9  "514) 

5 

FORMA  T 

(  1  A  * 

F  E  .  A 

,  2X  »  F 8 • 6 ) 

P  1 

FORMA  1 

(  )  HU 

*  F  X  * 

AE'Al  »  E  .6*  AH 

MIN.  ) 

P2 

FORMAT 

(  1  H 

*  sx  * 

F 1  c  .  ft  .  F  X  ,  ETF.ft, 

•  X  *  E  1  c  .  ft 

* 

f  X 

* 

p  P 

FORMAT 

(  1  H 

»  7  x  • 

A  F  1  0 . 2  ) 

P  U 

FORMA T 

(  1  H- 

*  TUX 

,  2  fl  HI >A  I  A  G  !  VF  N 

M F |  O W  ARE 

FOP 

0  S 

FORMA  I 

(  1  Hu 

*  F  X  » 

A  H 1  (  7  )  *  A  A  »  AMi  ( 

7 

)  .  F  X  *  AH 

T 

(  N- 

) 

1  AX  i  AMI 

(  N-2 

)  *  AX 

,  PHIS  ) 

P  A 

FORMA  1 

(  1  HU 

9  1  U  A 

,  /  H  1  G  *  !  H X  *  l  H X  * 

14X9  IMP  9 

1  EX 

• 

READ  ( S 

>  1  ) 

SK  * 

S»  i )  S  9  AR*  H  E  9  AS 
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V 

PF AD  ( F 

*  2  ) 

G  E  *  C 

P  1  9  CP/9  CP  B  9  El —  9 

Tr,  9  HR 

RFAD  (S 
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i  GO  * 

TO  9  PO*  XT;  9  D  1  *  F  E 

T 

)  *  H  e  o  (  i  7  ) 

F  P  E  (  2  A  )  9  Y  (  17) 
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F  F  A  0  (  A  *  4  )  I-  S  »  N  >  9  |\H  •  N4 

R  F  A  n  (5*5)  (  X  (  J  )  *  F  R  F  (  J  )  *  J  =  1  ,  M  2  ) 

WRITE  (5*54) 

WP  I  1  E  (  4  *  8  5  ) 

W  R  I  I  F  (5*86) 

AR=AR/1Z.0 

F M  8  =  N0- 1 

N  8  A  =  N  8  —  7 

N|8P=N8-4 

nR  =  A  R / FM8 

F  N  4  =  N  4 

N4  A  =  N4+8 

nZL=EL/FN4 

G  R  =  1 U • 7  8  * ]  .80 

1)1=01/54)0.0 

HR=HP*1  000.0 

F  I  F  =  5  0 . 0 

T I MF=0 . 00 

I K  =  S  K / ( 06*CS ) 

p I =8. 1 41 598 

A  A  1  =  8 . 0*6  K* (  I  , 0-V l/AW 

A  6  =  6 O R T  (  D  T  /  (  P  T*  I'  K  )  ) 

A8=-nT/AP+?,0*A5 

a  4  =  -n  T  * m  /  (  ? .  o *  a  r  )  4-4 .  n* n  i  *  a  4  /  8  .  o 
A  5  =  .S0R1  (  0  T  *  P  I  *  T  K  ) 
A7=I.j/(4.U*TK*nT)  ' 

A8=1.0/(2. U*SQR T ( I K*0 I ) ) 

B  3  =HF*AS/6F 

R?=0. 5*CPZ+3 . 5*CH i-( P] 

R  8  = V / ( GR*GF  ) 

R4= I  8 . 0 * D R * H F /SK 
F5=0. 000885 

D5= 1 ?.Q84 

R  7  =  1  Q  7  1  4 . 0 
88=4.858 
PO=SOQ?,U 
nO  9  T=i, M4A 

no  9  j= i , n 4 

9  T (  I , J ) =  TO 

00  10  t  =  l  ,  i\| Q 
F  I  =  T  - 1 

10  R (  T  )  =  F  T *0R 

00  11  I =1 , MR 

F'X  A  (  T  )  =0. 00 
FXr (  T  )  =0.00 

no  n  j= i *  nr 

P  A  (  I  «  I  )  =  0 . 0  0 
PF(  T,  |  )  =0.0  0 

11  PR (  I  *  J ) =0 . 00 
no  18  T  =1  ,N8 

U  A  =  R (  I  ) *R (  I  ) *A  7 

ur  = ( AR-R (  I  )  ) **Z*A7 

T  F  ( U A • G l  . T  IF)  nO  T O  17 
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♦  A  I 
*  .  •> ) 

•  ♦  q  '  . ,  f  =  a 
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.  s ) \  . r  =  r 

. 


’  *»/  •  r  e  T 

Of  c  f  L* T  )  I 

*•*  *  =  T 

.  rill 
fl**  f  =’ 

n  •  T  ) 

M  *  T  )  -  J 

f'1, 


1-10  9 


f-XA(  I  )=1.0/PXP(UA) 

IF  (  UC.GT  .  F  TF  )  GO  To  1  3 

p  xc  (  i )  =  i .  o/p  xp  ( ur ) 

1  3  TON  I  I  NUF 

no  ?0  T  =1  , N3 

no  ?n  j=i ,  i\n 

U A  = ( P (  I  )-P ( J )  )  **2* A  7 

UB= ( R (  I  ) +R ( J  )  )**?* A 7 

UF=  (  -2  •  0*AP+P  (  I  )  +  P  (  J  )  )  **2*  A  7 

IF  (  U  A  .  G  I  .  F  I  F  )  00  TO  2  2 

PA  (  I  *  J  )  =1  •  u  /  (-  XP  (  DA  ) 

2  2  IF  (UR.  or.  FTP)  OO  ro  2  4 

PR  (  I  , J ) =1  . u/ P  XP ( UP  ) 

2 4  IF  (Up.oI.PIP)  GO  ( 0  20 

PE (  I  * J ) =1 • 0/ FXP ( UF  ) 

20  fONTINUF 

no  110  T  =  1  ,  Nj  3 

X  J  = ( AP-R (  J  )  ) *AP 

119  DO  I  20  J= 1  ,N2 

IF  (XJ.L  F.X(J) )  GO  TO  121 
IF  (Xj.FQ.x(J) )  GO  10  122 
12  0  r  O  M  T  T  M  U  f 
XFP=0. 00 
go  ro  ii7 

12 2  XFR=1 . 0— FRF ( J ) 

GO  TO  117 

121  DX= ( XJ-X( J-l  )  )  /  6.0 
DO  123  JA=1,7 
F J A= JA-1 

X  X  = ( X ( J-l  ) +FJA*DX ) **2 
17  7  Y  (  J  A  )  =  1  .  ( j  /  F  X  P  (  X  X  ) 

X  f  p  =  Zt  i  .  0*  (  V  (  1  )  +  Y  (  7  )  )+2iP.  0*  (  Y  (  2  )  +  v  (  6  )  )  +2  ' . 0*  (  Y  (  3  )  +  v  (  e-  )  ) 
XFR=(XFP+?72.0*Y  (4)  )*nX*2.0/(  14  0. 0*.S  OP  I  (  P  T  )  ) 

XFR=1  .0-FPF(J-l  ) - X  F  R 
117  Rfc  R (  1  )  =XER 
1  1  6  TON  I  I NUF 

DO  73  1=1 ,N3 

F  X  A  (  I  )  =  R  (  1  )  *  P  (  1  )  *  F  X  A  (  f  ) 
no  73  J=1 ,N3 

73  H  A  (  I  ♦  J  )  =  (  P  A  (  I  *  J  )  -  P  R  l  1  *  J  )  -PP  (  UJ)  )  *  R  (  J  ) 

102  I  A  =  1 

T  T  MP  =  1  J  Mh +nl *X0 . C 
W R  1  T F  ( 6 , fl 1  )  I  IMF 

X!  =0.  10 
T  Y  (  1  )  =  T  OO 
T  Y  ( 2  )  =XO 
|  Y  (  3 ) =  pn 
1  0  3  P  =  1  (  T  A  ,  N  3  ) 

IF  (TA.Gi.4)  OO  10 

IP  (  1  A. to. 4 )  GO  T 0 
n/_  =nZL*0. 20 
GO  10  101 

■  73  P)Z  =DZL*0.  40 
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.  r*  ( i  *  I 

.  =  (  t  •  T 

,|  =  1  ^  f 

•  ' 

J  1  .  1  !  .  1  '  ) 

.  ■' 
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}  -  (  •  )  -*  M 

<  A  •  t '  .  A  1  ) 


i  -  j  in 


GO  10  30! 

76  I  )Z  =  i ) /  L 

1  o  1  A  3  =  A  A 1  *04 

I  b  =  R4*  r  Y  (  1  )  +  ]  R  .0*  |  (  I  A  (N'!-l  )  -V.  0*3  (  fA,|\H-2)+2.0*l  (  l  A  .N  4-3  ) 

Tb  =  F.S/  (  1  1  .  U+B4  ) 

0=  (  T  S-P ) / 01 
no  3  0  T  =  3  , M  3 

3  0  Y  (  T  )  =  R  (  T  )  *  I  (  1  A  «  I  )  *  (  —  H  P  P  (  T  )  ) 

W  3  =  u  .  u  *  Y  (  3  ) 

HO  Bui  k'  =  ‘S»M3R«4 
3  U  3  W  3  =W3+2.0*Y  {K  )+4*0*Y  (  K  +  2  ) 

A  2  3  =  I>P*2  .  0*  (  Y  (  3  )  +  Y  (  N 3  )  +  w  1  )  /  3  •  0 
W2=4.0*Y ( 2 ) 
no  B  0  2  K=B*N3A*2 
W2=W2+2.  )* Y ( K ) +4. 0*Y ( K  +  l  ) 

A  2  2  =  OR* ( Y (  ]  )+Y(NB)+w2)/3.0 
A2  6= ( A2  2- A2  3  )  /  3  3  •  0 
A2  = ( A22+A2  3 ) / ( AR*  I  K  ) 

r  b  =  (  (  3  .  0-  l  Y  (  2  )  )*rp3+u.c.*lY(2)*rpp  +  i.c;*TY(2)*rPB)*(TY(3  )  _  r  R  ) 

3  00  m  =  A  3  *  (  A  2  + A '3*P  +  A  4*n  ) 

F  (  3  )  =  T  Y  (  3  ) 

F  (  2  )  =  T  Y  (  2  ) 

F  (  B  )  =  j  Y  (  3  ) 

DO  40  K  =  I  »  4 

C  1  =  1  I  •  u -F  (2)  ) *F  XP ( R6-P  7/ F  (  1  )  ) 

(  2=u.b*F(2)*tXPlHH—  b  9  /  F  l  1  )  ) 

3  =  (  1  •  —  t  (  2  )  )  *  HI  +  •  b *F  (  2  )  *<  k/  +  j  .  -  *f  {  s  H'3 
S(2*K)=R3*l)2*F(3)*(Cl-C2)*B600.0/(t-M  )  *  (  i  .  U  +  F  (  2  )  )  ) 

S(  3  9  K  )  st nZ *  (  R 1  *  {  H+Q*0  «  f*d T ~F  (  i  >  )  —  (  B2*  (  H  (  1  )  —  T B  )  +HR  )  * S  (  2  *  K  )  )  /  (  3 
S  (  B,<  )  =-DZ*{  3  .  0  +  F  (  2  )  )  *  I-  (  I  )*R2/l-  l  B  ) 

IK  ( K . G T . ? )  GO  3  0  4  2 
F (  1  ) = I Y (  I  ) +  b (  1  9<  )  / 2 . 0 
M2)  =  IY(2)+M2»K)/2.U 
F  (  3  )  =  I  Y  (  3  1  +  b  (  3  *  K  1/2.0 
GO  10  40 

A  2  ft-  (  K  .G 1  .  3  )  GO  I  O  44 

F(1)=IY(1)+S(1*K) 

F!2)=lY(2)+b(2»K) 

F  (  3  )  =  I  Y  (  3  )  +  b  (  B  ,  K  ) 

GO  10  40 
40  TONI fNUF 
44  no  4b  T  =  3  , B 

4b  F  (  T  )  =  f  Y  (  II  +  (  S  (  I  »  3  1  +b  (  T  9  U  )  +2,  •  0*  (  S  (  ♦  2  ) +b (  1  . 3 )  )  ) / 6 . 0 

04=  (  (  3 . 0-F  (  2  )  )*rP3-t-0.^*M2)*rp?  +  3.£*c-(2)*rp'-M*(M3  )-  F) 
n2=-oc*n | * ( r4-rc+ ( F ( 2 1 - 1 Y ( 2 ) ) *HP l 
f  F5  =  APS ( O  3 -0  2 1 
AQ= Abb ( 01 +02 ) 

I  F  b= IF  6*2 .  0  /  AO 
IF  (fF_S.LF.FS)  GO  TO  bo 
0=1  (03+02*2.0)/  l  3 • O* A  J  )  - A  2 - A  3  *  P )  / A  4 
GO  i O  3  Ou 
bo  no  b  6  T=2.N3 
no  b  6  J  =  3  .MB 
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A*  fat 

•  f  )  =  t 

- 

a  (  »  f  ) 
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•  H  M  l 

' 


56  V  ( j )  =  i  (  IA  » J ) *RA (  I 9 J ) 
iA,  i  =  4 .  u*y  (  3  ) 
no  3  9  K  =  6»iM^k,4 

W 1  =  W I 4-  2 • U *  Y  ( K  ) +4 • 0*  Y l *  + /  ) 

A  ?  1  =  OP* ?  .  0*  (  Y  (  1  )  +  Y  (  N 3  )  +  v4 1  )  /  3  •  0 
W  ?  =  4 . 0  *  v  (  7  ) 
no  6 o  |Y  =  3.N7A.? 

9  Q  W?  =  W?4-2  .  0*  Y  (  <  )  +4 . 0*Y  (  v  +  1  ) 

A  7  7  =  no* ( y  (  1  ) + Y ( N 9 ) +  W?  )  /  3  . 0 
A?3=(A??-A?1  )/] 5.u 
T  A  =  (  A  ?  /  +  A  2  3  )  /  (  /  •  u  *  R  i  1  )  *A  5  ) 

T  A  1 =  AR*R*BF_R i  |  )  / K (  {  ) 

TR=(DT+(AR-P(  |  )  )**//  t  TK  +  TK )  )  *-R F  K  (  i  ) 

T  C  =  -  (  A  R  -R  (  i  )  )  *  A  6  *  F  X  (  (  T) 

TCI  =  ( TR  +  TC ) *AR*Q/R (  1  ) 

6  A  T  7  (  I  )  =  1  A  +  T  A 1  +  T  r  1 
TT ( 1 ) =1 600. r 
T  Y  (  l  )  =  F  (  1  ) 

T  Y ( ? ) =F ( 2 ) 

T  Y  (  9  )  =  F  (  9  ) 

nb  T  =o*n7  +  P 

XL-XL+OZ 

W P  I  l  F  (  6  f  P  9  )  I  i  (  9  )  9  !  I  (  7  )  »  T  T  (  N  3 - 6  )  *  ]  i  (  i\j  3 - 4  )  ,  T 7  (  N  9 -  2  )  »  l)F 
W  R  I  I  F  (  6  »  P  2  )  1  Y  (  1  )  »  1  Y  (  2  )  «•  7  V  (  3  )  »  XL 
IF  (IA.GF.fM4A)  GO  TO  61 

HQ  6  4  T  =1  ,  N3 

64  T  {  I  A  *  I  ) =T T (  I  ) 

I A= I A+i 
GO  TO  1 09 

6i  tf  (TTMF.|_T.FrT  )  no  70  1  07 
roMT T  MUF 
FMD 

$  F  M  7  P  Y  FIMA!. 

4.R60  0.403?  ?  3 1  .  8  P  0.50  8  3.00  8  6.^0  0.40  0 

)•  00  .  .  4.  .  1673.0  26. 

Q4o.UF  ?  J  0  0 . 0  0  16.70  7.00  5.00  3. OQ 
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PART  TWO 


CORRELATION  OF  MULTIVARIATE  DATA 
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I .  INTRODUCTION 


The  task  of  finding  suitable  approximating  functions 
may  arise  both  in  data  correlation  and  in  data  representation. 
It  is  assumed  that  the  approximating  function  is  a  linear 
combination  of  a  set  of  specified  functions  and  that  the 
coefficients  are  to  be  determined  using  some  criteria  of 
best  fit. 

For  data  correlation,  the  form  of  the  functions  may 
be  known  but  the  data  points  which  are  to  be  fitted  may  con¬ 
tain  errors,  where  these  errors  are  randomly  distributed. 

Then  the  criterion  of  least  squares  is  normally  used. 

For  the  data  representation  the  data  points  are 
assumed  to  foe  free  of  error  and  our  task  is  to  find  an  analy¬ 
tical  approximating  function,  which  may  be  more  convenient  for 
computational  purposes.  In  this  case,  the  Chebyshev  criterion 
may  be  more  appropriate  as  a  criterion  of  best  fit. 

In  general  the  problem  may  be  stated  as  follows: 

It  is  desired  to  approximate  a  real-valued  function 
u.  by  n+1  known  functions  of  g.(x.  ,  y .  ,  z .  ,  ...)  for  0<j<n, 

2.  Jill 

at  a  set  of  points  (xj ,  yj ,  z^,  ...)  for  i  =  1,  2,  ...,  m, 

continuous  on  a  closed  region  R.  If  the  deviation  at  each 
point  is  defined  as 


n 

^  *  Aj  g  j  (xj,  ^  j '  •••)  ( I  •  1 ) 

j-0 


where  x,  y,  z  €  R 
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then  Aj  are  the  n+1  expansion  coefficients  to  be  evaluated 
according  to  the  following  criteria: 

1.  least  squares  criterion 


m 

V  2 

i.e.  )  is  a  minimum 

i=l 


and 


2 .  Chebyshev  criterion 


i.e. 


€  . 

i  max 
i 


is  a  minimum, 


These  are  by  no  means  the  only  criteria  for  approxi¬ 
mation.  Indeed  there  has  been  considerable  work  done  in  develop 
ing  the  "Least  Qth"  method(l),  for  which  the  "least  squares" 
is  but  a  special  case.  If  we  define  the  function  S  by 


m 


i=l 


(q  integer) 


(1.2) 


and  it  has  only  one  real  minimum,  then  it  is  required  to  find 
the  expansion  coefficients  Aj  in  minimizing  S.  When  q  is  any 
integer  other  than  one,  the  system  becomes  implicit  in  Aj  so 
that  the  method  of  steepest  ascent  or  other  trial-and  error 
approaches  must  be  used  in  solving  the  problem.  The  com¬ 
plexity  and  the  often  slow  rate  of  convergence  are  the  draw¬ 
backs  in  rendering  the  criterion  of  "least  Qth"  rather  unpopular 
When  q  is  sufficiently  large  it  can  be  seen  that 
the  criterion  of  "least  Qth"  approaches  that  of  Chebyshev 
approximation  in  minimizing  the  maximum  deviation. 
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II.  CURVE  FITTING  BY  LEAST  SQUARES  CRITERION 


The  least  squares  principle,  originally  proposed 
by  Gauss,  postulates  that  the  best  fit  occurs  when 


m 


i=l 


(II. 1) 


is  a  minimum,  where  w^  is  the  weighting  factor  associated 
with  each  data  point.  For  a  function  of  an  independent 
variable  x,  to  be  approximated  by  some  general  polynomials 
Gj(x),  at  most  of  degree  j,  equation  (II. 1)  may  be  written 
as 


S 


m  n 


i=l  j=0 


(xi) 


(II. 2) 


Since  S  is  a  quadratic  function  of  A j ,  its  minimum  exists  at 


dS 

-  =  0  for  j  =  0,  1,...  n  (II. 3) 

SA3 

After  expansion  it  yields 

n  m 

y  a.  y  w.-g.  (xi)-ck  (x.) 

j=0  i=l 


m 


i  =  l 


wi  ui  °k  (xi> 


k  =  0,  1,  2, 


(II. 4) 


n  . 
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If  we  let 


m 


wi  Gj  (xj_)  Gk  (xi) 


(II. 5) 


i=l 


and 


m 


z 


i=l 


then  equation  (II. 4)  can  be  written  algebraically  as 


DA 


E 


(II. 6) 


This  is  a  system  of  n+1  linear  equations.  When  the  matrix  D 
is  non-singular  it  follows  that  Aj  is  unique  for  minimizing 
S.  The  evaluation  of  Aj  in  equation  (II. 6)  can  be  achieved 
by  any  standard  method  in  the  field  of  linear  algebra. 

To  check  the  closeness  of  fit,  the  criterion  of 
variance  8  may  be  used,  which  is  defined  as 


2 


i=l 


(II. 7) 


5 


m  -  n 

The  variance  is  a  measure  of  the  average  deviation  of  the 
experimental  points  from  the  approximating  function.  The 

t 

is 

degree  of  the  polynomial  n a  arbitrarily  chosen  at  some  approp¬ 
riate  level,  however,  it  must  be  less  than  m.  In  going  through 
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the  analysis,  the  variance  is  calculated.  Then  n  is  increased 
by  one  and  the  analysis  is  repeated,  until  there  is  no  signifi 
cant  change  in  the  variance.  However,  a  complete  set  of  cal¬ 
culations  has  to  be  carried  out  for  each  additional  polynomial 
tested;  and  the  Sets  of  expansion  coefficients  determined  bear 
no  relationship  to  the  others.  Furthermore,  computational 
difficulties  arise  from  the  fact  that  the  matrix  D,  in  most 
cases,  is  ill-conditioned  -  because  of  the  sensitivity  of  Aj 
to  slight  change  of  values  in  D, 

The  solution  to  both  these  problems  -  the  ill 
conditioned  matrix  and  the  necessity  for  recalculating  the 
coefficients  for  each  higher  degree  polynomial  tested  -  lies 
in  the  use  of  orthogonal  polynomials. 

When  a  sequence  of  polynomials  J  (x)  j  is  ortho¬ 
gonal  in  the  summation  sense  over  a  point  set,  they  satisfy 
the  following, 

m 

w- Gj  (xi)-  G-k  (x^  =  0  j  /  k 

i=l 

m 

=  w- Gj  (xj)  j  =  k.  (II. 8) 

i*=l 

When  this  is  substituted  into  equation  (II. 4)  the  expansion 
coefficients  can  be  evaluated  very  easily  as 
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A 


m 

X  wiui'Gj  (x^ 

i=l 


w 


(xi) 


(II. 9) 


without  the  necessity  of  matrix  manipulation.  Furthermore, 

the  coefficients  are  independent  of  one  another  and  of  the  value  n. 

With  the  orthogonality  of  (II. 8)  and  the  expression 
(II. 9),  equation  (II. 2)  yields 


m  n 

S  =  >  w.  J  u2.  -  2u.  V  A .  G .  ( x .  ) 

L  1  ]  1  1  L  3  ]'  i 

i=l  1  j=0 


+ 


h  \ 

Aj  Gj  (Xi)  ^  Ak  Gk  (Xi)  | 

k=0 


m 


i=l 


n  m 

j  X  Wi  Ui’  G3  (Xi> 

j=0  i=l 


n  n 


m 


+  ^  \  y  wi'Gj  (xi)  Gk  (xi) 

j=0  k=0 


i=l 


m 


i=l 


2 


n 


j=0 


i=l 
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n 


m 


+ 


£  wA  Gj2  (x±) 
i=l 


m  n  m 


i=l  j=0  i=l 


(11,10) 


Since  S,  by  definition,  is  positive  definite,  as  n  increases 
S  must  decrease.  This  verifies  the  convergence  of  the 
approximation . 

The  literature  on  the  subject  of  orthogonal  poly¬ 
nomials  is  voluminous.  However,  we  need  be  concerned  with  only 
certain  properties.  Thus,  all  orthogonal  polynomials  satisfy 
the  recurrence  formula(2) . 


Pn+l(x) 


(A  x  +  B  )  P  (x)  -  C  P  _  (x) 
n  n  n  n  n-1 


n  =  1,  2,  3  ...  (II. 11) 

where  An,  Bn,  Cn  are  arbitrarily  constants  to  be  constructed. 
In  the  following  calculation  AR  is  taken  to  be  unity  as  sug¬ 
gested  by  Forsythe(3),  so  that  the  complete  recurrence  form 
is  expressed  as  follows: 


» 
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PQ(x)  =  1 

Px(x)  =  (x  -  ax)  Pq  (x)  (11.12) 

P2(x)  =  (x  -  a2)  Px  (x)  -  Pq  (x) 


P  j  (x) 


(X  -  a.)  P^Cx)  -  PjPj_2(x) 


where  a_., 


a  . 
D 


are  defined  as  follows: 
m 


i=l 


m 


i=l 


(x±) 


(II. 13a) 


m 

X  wi  xipj-l 

i-1 


(xi>  Pj_2  (Xi) 


m 


i=l 


(II. 13b) 


They  can  therefore  be  calculated  sequentially  in  a  systematic 
manner.  Their  derivation  and  proof  of  orthogonality  are  given 
in  Appendix  A. 

An  advantage  of  constructing  such  polynomial  over 
other  orthogonal  polynomials  such  as  Chebyshev,  Gram,  Laguerre 
and  many  others,  is  that  it  does  not  require  specified  spacing 
of  the  data  points. 
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A  study  of  the  effect  of  weighting  factors  has  also 
been  investigated.  By  modifying  the  weighting  factors  linearly 
according  to  the  absolute  deviations  obtained  in  the  analysis 
using  equal  weighting  factors,  it  was  found  that  the  maximum 
deviation,  though  not  necessarily  the  variance,  can  be 
reduced . 

Numerical  illustrations  are  given  in  a  later  section. 
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III.  SURFACE  FITTING  BY  LEAST  SQUARES  CRITERION 


This  is  an  extension  of  curve  fitting  in  section  II. 
Here  we  try  to  find  the  approximating  function  for  systems 
with  two  independent  variables.  The  approximating  function 
to  be  investigated  is  a  mixed  polynomial  of  the  forms 


n  m 

f(x,y)  =  y  y  Aij-Pi  (x)  •  Q j (y)  (III.l) 

i=0  j=0 


where  n,  m  are  the  highest  degree  of  the  polynomials  in  x 
and  y,  respectively. 

P(x) ,  Q(y)  are  the  orthogonal  polynomials  to  be 
generated  from  given  data  through  recurrence  formulas  as 
shown  below: 

At  y  =  constant 


and  at  x 


po(x) 

=  i 

PX(x) 

=  (x  - 

-  0^) 

PQ(x) 

P2(x) 

=  (X  - 

■  a2] 

Px(x)  - 

•  •  • 

P^x) 

=  (x  - 

-  a±> 

Pi-l(x) 

constant 

Q0<y) 

=  1 

Qi(y) 

=  (y  - 

-  7i> 

Q0(y) 

q2(y) 

=  (y  - 

'  y2> 

Q^y)  - 

Qj(y) 

=  (y  ■ 

-  Yj) 

Qj-i(y) 

e2  P  (x)  (III. 2) 

-  Pi-2(x> 


62  Qo(y) 

-  6j  Qj_2(y). 


(III.  3) 
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If  the  data  point  at  (x^/Y^)  is  z]^  i*  then  the  error 
of  approximation  is  given  by 


£  -i  i  =  z, 
k,  1  k 


,i  -  f(vyi>. 


(III. 4) 


For  a  rectangular  network  with  u  and  v  numbers  of  points  along 
the  x  and  y  axes  respectively,  the  sum  of  the  error  squares 
becomes 


S  = 


u 


k=l  1=1 


(  £  )2 

k,l 


u 


V 


n 


m 


=  I  I  i  zki  - 1  I  Aijpj(xk)Qj(yi)}2 

i=0  j=0 


(III. 5) 

2 


k=l  1=1 


Here  weighting  factors  are  not  considered  and  S  is  a  function 
only  of  A^j.  By  the  criterion  of  least  squares,  it  is 
necessary  to  find  the  matrix  (A^j)  such  that 


dS 


dA 


st 


=  0  for  s  =  0,  1,  . . . ,  n 

t  =  0,  1,  . . . ,  m 


(III. 6) 


After  expanding  equation  (III. 6)  and  rearranging, 


this  gives 


2-12 


n  m  u  / 


— i  •;  n 


Zi  L  Zj 

i=0  j=0  k=l  1=1 


Pi(xk>  Ps(xk)Q.(yi)  Qt(Yl) 


u  v 


X  Yj  Zkl  PS(Xk)  Qt(yi} 


(III. 7) 


k=l  1=1 


Now  if  P^(x)  Q.(y)  are  orthogonal  polynomials  in  the 


summation  sense,  i.e. 


u  v 


X  X  pi(*k)  Ps(xk)  Qj<yi>  Qt(n> 


k=l  1  =1 


=  0  if  i  ^  s  or  j  t/-  t 


but 


u  v 


X  Ps2(xk>  Qt2(yi> 


k=l  1=1 


if  and  only  if 
s  =  i  and 

t  =  j  (III. 8) 


then  the  expansion  coefficients  can  readily  be  evaluated  as 


u  v 


Y  X  Zki  PiUk)  Qi<Yi) 


k=l  1=1 


A.  .  = 
ID 


(III. 9) 


u  v 


X  I  Pi2(xk)  Qj2(n> 


k=l  1=1 
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and 


It  has  been  proved  in  Appendix  B  that  if 


u 


^  xk  pi-l(xk) 


k=l 


a . 
1 


P.  = 

1 


u 


I  pLi(v 


k=l 


u 


Z  **  Pi-l^xk)  Pi-2(xk) 

k=l 


u 


Z  Pi-2(xk> 


k=l 


V 

z  yl  Qj-l(yl> 


1=1 


U 


Z  °j-i(yi: 


1=1 


(ill. 10a) 


(III. 10b) 


(III. 10c) 


Z  yl  Qj-I(yi} 


1=1 


Z  °j-2(yi) 


1*1 


(Ill.lOd) 
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then  Qj(y^)  will  indeed  satisfy  the  orthogonality 

condition  of  (III. 8) 

The  application  of  this  analysis  has  given  very 
satisfactory  results  in  comparison  with  that  using  Chebyshev 
polynomials.  However,  there  is  difficulty  in  incorporating 
any  weighting  factors.  Furthermore,  this  approach  requires 
the  data  points  to  be  on  a  rectangular  grid. 
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IV.  FUNCTION  APPROXIMATION  BY  CHEBYSHEV  CRITERION 


The  statement  of  the  problem  is  as  follows:  Given 
a  set  of  real  functional  values  f(x)  at  x  =  x^,  (i  =  1,2, ...  m) 
and  n+1  known  functions  gj(x)  for  0<j^n,  continuous  on  a 
closed  region  R,  find  an  n+1  dimensional  real  vector  a  =  (ctj) 
such  that  the  norm  of 


A  =  f  ( x ) 


Oj  gj(x) 

x  £  R 


( IV  .  1 ) 


is  a  minimum. 

This  problem  can  be  solved  by  some  iterative  method (4), 
based  upon  the  theorem  of  de  la  Vallie  Poussin(5),  which 
states  that 


(i)  The  solution  vector  a j  which  furnishes  the 
best  approximation  to  (IV.l)  gives  at  least  n+1  residuals 
equal  to  IM, 


where 


M 


minimum 
a  _i 


n 


f ( *i )  -  )  ctj  gj(x±) 


j=o 


max 
i 


(ii)  The  minimal-maximum  residual  M  of  a  system 
of  m  linear  equations  in  n+1  unknowns  is  equal  to  the 
largest  of  the  minimal-maximum  residuals  Mj_,  M2,  ...  cor¬ 

responding  to  the  various  systems  which  can  be  found  in 
taking  n+2  equations  at  a  time. 


*  r 


2-16 


In  general  the  method  can  be  described  as  follows: 

1.  choose  arbitrarily  n+2  points  f(x^)  with  R(  forming  a 

set  of  n+2  simultaneous  equations, 

2.  solve  the  simultaneous  equations  to  yield  aj"^  and  A^^ 
(where  the  superscripts  represent  the  number  of  trials) , 

3.  find  a  point  f(xr)  in  R  other  than  those  in  the  set 


f(x  ) 
r 


a 


g .  (x  ) 
3  r 


>  * 


(1) 


4.  replace  one  of  the  points  f(x^)  by  f(xr)  to  form  a  new 
set, 

5.  now  repeat  the  calculations  until 


f(xi) 


n 


j=0 


for  i 


1,2, 


m 


At  worst  this  method  requires  a  selection  and  solution 
of  (n+2)  systems  of  n+2  equations  in  n+2  unknowns.  When  m  is 
large,  which  is  usually  the  case,  the  method  becomes  so  time 
consuming  as  to  be  impractical . 

Since  the  problem  is  to  find  the  minimum  A  subject 

to 


n 

f(xjL)  -  y  cij  gj(xi) 
j=0 


<  A 


1  ^  i  <;  m 


d-m 


2-17 


i  .e. 


n 

X  +  Yj  aj  gj^Xi*  ^  f(^i) 

j=0 


X 


n 

Y  a j9 j (xi)  <  -f(xi) 

j=0 


(IV. 2) 


it  can  be  framed  as  a  linear  programming  problem. 
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V.  LINEAR  PROGRAMMING  -  SIMPLEX  METHOD 


It  is  appropriate  here  to  give  a  brief  description 
of  linear  programming  and  its  solution  by  the  Simplex  method. 
Detailed  discussion  of  the  method  can  be  found  in  many  standard 
textbooks (6, 7,  8) . 

The  general  linear  programming  problem  is  to  find 
a  vector  a  =  a(a-^#  ,  an)  which  minimizes  the  linear 

form  (i.e.  the  objective  function) 

n 

c  j  a  j  =  B  (V .  la) 

j=l 

subject  to  the  linear  constraints 


a  . 
D 


0 


^  *•  1 1  2 1  ♦  «  «  n 


(V.lb) 


and 


n 


3~- 


a.  .  ot  •  ~  b  ■  i  —  1  /  2  f  ...  m 

ID  3  i 


(V.lc) 


where  the  a. b.  and  c.  are  given  constants  and  n  <  m. 
ID  1  D 

If  the  constraints  involve  inequalities,  slack 


variables  s^  are  introduced. 
Therefore, 


ailal  +  ai2a2  + 
may  be  written  as 

ailal  +  ai2a2  + 
where  s^  is  non-negative. 


+ 


ainan 


•  •  • 
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If  some  variables  are  not  constrained  to  be  non- 
negative  then  it  is  necessary  to  let 


a 


1 


a 


(2) 

j 


because  any  number  can  always  be  expressed  as  the  difference 
of  two  non-negative  numbers. 

If  the  objective  function  is  required  to  be  maximized 


n 


j=l 


it  can  be  rewritten  as 


max . 


j=l 

By  arithmetic  manipulation  the  given  system  (IV.1) 
can  readily  be  transformed  into  a  canonical  form  as 


al 

al,  m+lam+l 

+ 

...  + 

a '  a 
l#n  n 

r-H 

il 

a2 

a2 , m+lam+l 

+ 

...  + 

al  a 
2,n  n 

-  (N 
12 

II 

•  •  « 

(V.2) 

a  + 

m 

am, m+lam+l 

+ 

...  + 

a 1  a 
m,  n  n 

II 

cm+l  am+l 

+ 

...  + 

c '  a_ 

n  n 

=  B  - 

Bo 

where  BQ  is  a  constant. 

Here  a1#  a2,  •••  #  am  are  the  basic  variables. 

Hence,  the  canonical  form  for  any  basis  is  the  rearrangement 
of  the  original  system  so  that  the  variables  become  unit  vectors. 
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The  Simplex  method  is  an  iterative  procedure  which 
searches  for  the  optimum  from  one  basic  feasible  solution  to 
another  basic  feasible  solution,  having  a  lower  value  of  B* 
Thus,  the  Simplex  method  deals  only  with  the  canonical  form  of 
bases  and  consists  of  the  following  steps: 

1.  examine  the  c and  find  the  smallest  one  (say)  c', 

D  s 

2.  if  c'>0,  then  the  current  solution  is  optimal, 

s 

3*  if  c'<0,  from  the  column  of  coefficients  a!  pick  out 

o  JL  o 

those  a!  >0, 
is 

4.  if  all  a!  s*  0,  then  no  unique  optimal  solution  exists, 

3.  s 

5.  if  some  a!  >  0,  find  the  i  for  which  b! /a!  is  the  smallest 

xs  1  is 

and  let  this  i  be  r, 

6.  using  a|g  as  the  pivot  element,  obtain  the  new  canonical 
form  in  which  a  has  replaced  a  as  a  unit  vector. 

S  j- 

7.  Repeat  the  whole  process  until  all  c^  >  0. 

In  many  cases  an  immediate  basic  feasible  solution 

cannot  be  obtained,  so  artificial  variables  (an+^,  an+2'***' 

a  ,  )  must  be  introduced.  Their  removal  is  accomplished  not 
n+k 

by  considering  at  first  the  original  objective  function,  but 
by  minimizing  the  infeasibility  form  W,  defined  by 

W  =  M  (an+1  +  an+2  +  ...  +  an+k)  (V.3) 

where  M  is  a  large  penalty  coefficient. 

The  minimization  of  (  V.3)  is  called  "phase  I", 
while  the  minimization  of  the  original  objective  function  is 
called  "phase  II" .  The  process  from  one  feasible  solution  to 


< 
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another  is  sketched  schematically  in  Figure  1. 

The  coefficient  dj  can  be  regarded  as  similar  to 
that  denoted  by  c^  and  associated  with  the  infeasibility  objec¬ 
tive  function  of  (V  .3) .  When  all  d^  ^  0  and  W  =  0,  end  of 
phase  I  is  reached.  But  if  W,  at  minimum,  does  not  equal  to 
zero,  no  feasible  solution  exists. 

The  overall  procedure  in  solving  linear  programming 
by  Simplex  method  is  illustrated  in  the  flow  diagram  of  Figure  2. 

In  function  approximation  the  problem  is  to  find  the 

minimum  X 


l 


.e . 


f  (x^ 


a. 


gj(xi) 


for  1  ^  i  <  m 


X  + 


n 

^  o u  gj(x±)  >  f(x±) 

j=0 


With  the  introduction  of  slack  and  artificial 


variables  the  problem  becomes 
to  minimize 

Z  =  X  +  M 


m 


i=l 


u . 

l 


(  <  (j  7  + 


<  J  :  - 


k+1  iteration  k  iteration 
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Figure  1 

Diagrammatical  Representation  of  Simplex  Method 
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Flaw  Diagram  for  Linear  Programming  by 

the  Simplex  Method. 


Figure  2 
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(continued) 
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subject  to 


n 


+  X  (ylj  -  y2j>  gj{xi>  -  vli  +  ui 
j=0 


f  (xi) 


*  - 


n 


X  (ylj  '  y2j>  gj(xi>  +  v2i  =  -f<xi) 


j=0 


1  <  i  <  m 


(V.4) 


where  v^,  w2i  are  t^le  s^ac^  variables  and  are  the  artificial 
variables  with  penalty  coefficients  M,  However,  the  variables 
cij  are  not  restricted  in  sign,  so  put 

aj  =  *lj  *  *2j. 

Thus,  all  the  variables  y^,  y2j#  v2i'  vi  and  ui  are 
non-negative. 

It  is  obvious,  therefore,  the  new  problem  may  be 
too  big  for  computation  if  m,  the  number  of  points,  is  large. 

An  elegant  way  of  getting  around  such  a  difficulty  is  pro¬ 
posed  in  the  following  section  by  applying  the  dual  linear 
programming  technique. 
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VI.  FUNCTION  APPROXIMATION  BY  DUAL  LINEAR  PROGRAMMING 

By  definition  if  two  linear  programming  problems 
are  said  to  be  dual,  they  are  related  to  one  another  as 
follows : 

when  the  minimizing  problem  is  considered  as 
primal  such  that: 

minimize  B  =  c,x,  +  c~x~  +  ...  +  c  x 

11  2  2  n  n 

subject  to 

n 

V  a.  .x.  >  b.  ( 1  <  i  <.  m)  (VI. 1) 

j=l 

with  all  x  j  >  0  ( 1  <  j  <  n) 

then  maximizing  problem  will  be  considered  as  dual  such 
that 

maximize  Z  =  b^y-^  +  +  •  *  *  +  kmym 

subject  to 

m 

^  aji^i  <  cj  (l<j<n)  (VI. 2) 

i=l 

with  all  y 0  (l<i<m) 

The  relationships  between  the  two  systems  are  such 

that: 

(i)  if  the  objective  function  of  one  is  to  be  minimized,  then 
that  of  the  other  is  to  be  maximized, 

(ii)  the  coefficient  matrix  of  one  system  is  the  transpose  of 


the  other . 
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(iii)  The  constants  on  the  right-hand  side  of  one  system  are 

the  coefficients  of  the  objective  function  of  the  other. 

As  a  general  case  consider  a  primal  to  be  a  mixed 
system  consists  of  the  following: 

(a)  k  constraints  of  type  ^  aijxi  =  bj  j  =  1,  k 

(b)  r  constraints  of  type  ^  aijxi  <  bj  j  =  k+1, ...,k+r 

(c)  s  constraints  of  type  a^x.^  >  bj  j  =  k+r+1, . . . ,k+r+s 

Some  of  the  variables  are  non-negative  while  the  rest  are  un¬ 
restricted  in  sign.  This  primal  system  is  the  most  general 
system  possible.  It  is  represented  diagrammatically  in 
Figure  3,  with  letters  A,  B,  C,  D ,  E  and  F  denoting  the  sub¬ 
arrays  of  the  coefficient  matrix. 

Before  dualizing  the  primal,  it  is  necessary  to 
multiply  the  constraints  k+1,  k+2,  ...,  k+r  by  -1  so  as  to 

reverse  the  inequalities.  The  dual  is  shown  schematically  in 
the  lower  diagram  of  Figure  3,  with  superscript  T  denoting  the 
transpose  of  the  original  subarray  of  the  coefficient  matrix. 

The  fundamental  property  which  is  relevant  in 
function  approximation  can  be  stated  as  follows: 

If  X  is  a  feasible  solution  vector  to  (VI. 1)  and 
Y  a  feasible  solution  vector  to  (VI. 2) ,  and 

n  m 

bi  Yi  (VI. 3) 


emc 


Figure  3 

Duality  of  Linear  P rogramming  Problems 
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then  X  is  an  optimal  solution  vector  to  (VI.  1)  and  Y  is  an 
optimal  solution  to  (VI. 2)  . 

Now  if  we  consider  the  system  (IV. 2)  again,  since 
are  not  restricted  in  sign,  the  dual  system  with  reference  to 
Figure  3,  may  be  written  as  follows: 

ff(ui  -  vi)  =  Z  (VI. 4a) 

(ui  +  Vi)  <  1  (VI. 4b) 

gij'^ui  “  =  0  (VI. 4c) 

and  j  =  1,2,  ...,  n. 

Now  to  convert  the  constraint  (VI.4b)  into  an 
equation,  a  slack  variable  has  to  be  introduced.  However, 
if  we  assume  that  we  are  dealing  with  the  non-trivial  case, 
in  which  A>0,  this  slack  variable  will  not  enter  the  optimal 
solution  at  a  positive  level (14)  -  in  other  words  this  slack 
variable  cannot  be  in  the  final  basis  and  must  therefore  take 
a  value  of  zero.  Thus  the  constraint  (VI. 4b)  can  be  written 
as 


subject  to 


and 


m 


i=l 


m 


i=l 


m 


i=l 


for  all  Ui ,  Vi >  0 


(u .  + 


vi) 


i=l 


=  1 
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When  the  optimal  solution  is  reached 


(B) 


(A) 


opt .  '  ' ' opt . 

it  can  be  seen  that  the  optimal  basic  vector  of  the  dual 
system  will  correspond  to  those  constraints  in  the  primal 
that  gives  the  optimum  A .  After  solving  these  constraints 
in  the  form  of  equalities  the  various  coefficients  ctj  and  A 
will  be  determined. 

The  advantage  of  this  approach  is  that  the  number 
of  constraints  in  the  dual  is  equal  to  the  number  of  coef¬ 
ficients  plus  A;  thus,  it  is  independent  of  the  number  of 
data  points.  Furthermore,  if  an  extra  coefficient  in  the 
primal  system  is  required,  it  is  equivalent  to  an  additional 
constraint  in  the  dual. 

As  stated  in  equation  (1.1),  g. .  can  be  any  known 
functions  with  any  number  of  independent  variables.  More¬ 
over,  the  data  points  can  be  located  in  any  convenient  or 
random  spacing,  which  is,  indeed,  a  very  significant  advan¬ 
tage  over  other  methods  of  approximation  involving  more  than 
one  independent  variable.  This  approach,  in  fact,  enjoys  a 
high  degree  of  freedom  in  approximation. 


VII.  APPLICATION  OF  DUAL  LINEAR  PROGRAMMING 


TO  MULTIPLE  REGRESSION  ANALYSIS 

In  order  to  obtain  a  mathematical  description  of 
the  relationship  between  the  functional  value  and  its  variable 
there  are  several  ways  of  approaching  the  problem, depending 
on  the  characteristics  of  the  observations  and  on  the  aim  of 
the  analysis.  Correlation  analysis  and  regression  analysis 
are  two  of  the  statistical  methods.  The  former  has  been 
discussed  in  the  previous  section  while  the  latter  is  to  be 
described  briefly  in  the  following. 

The  purpose  of  regression  analysis  is  to  obtain  a 
best  fit  of  a  set  of  observations  of  independent  and  dependent 
variables  by  an  equation  of  the  form 

n 


y 


i=0 


where  y  is  the  dependent  variable,  x^,  X2 ,  ...  are  the 

independent  variables  in  some  given  functions, fi  and  bi 
are  the  coefficient  to  be  determined.  However  it  is  not 
known  whether  all  the  independent  variables  x^  and  the  given 
functions  f^  will  provide  any  significant  contribution  to  the 
approximation.  It  is  the  purpose  of  the  analysis  to  test  tneir 
significance.  The  criterion  for  accepting  fi  to  be  an  approxi- 
mant  is  the  extent  of  improvement  in  the  variance  upon  the  in¬ 
clusion  of  f^  into  the  approximating  equation.  Hence,  in  the 
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course  of  analysis  if  an  independent  variable  x^  or  a  given 
function  fj_  is  tested,  a  least  square  analysis  is  performed  and 
the  variance  is  calculated.  Obviously  for  multiple  regression 
this  method  is  very  time  consuming. 

A  stepwise  procedure  has  been  proposed(9),  which 
introduces  a  variable  or  a  given  function  at  a  time.  However, 
a  least  squares  analysis  must  be  used  in  obtaining  the  coefficients 
before  the  variance  can  be  determined. 

If  the  Chebyshev  criterion  and  not  the  variance  is 
to  be  used  as  a  criterion  for  rejecting  the  variables  or  given 
functions  and  for  the  goodness  of  fit,  then  dual  linear  pro¬ 
gramming  can  be  used  in  the  analysis.  Here  we  determine  whether 
or  not  there  is  a  significant  change  in  the  maximum  absolute 
deviation . 

As  seen  from  the  preceding  section,  with  an  intro¬ 
duction  of  an  approximating  function  fj_,  an  extra  expansion 
coefficient,  b^  is  to  be  determined,  which  is  equivalent  to  an 
addition  of  a  constraint  in  the  dual  problem.  Corresponding 
to  each  additionafof  constraint  to  the  system,  an  optimum  ob¬ 
jective  function  (B)  t  can  obtained  which  is  the  maximum 
absolute  deviation  of  the  system.  Hence,  to  see  the  significance 
of  a  certain  approximating  function  f^  it  is  necessary  only  to 
compare  the  values  of  the  objective  function  before  and  after 
the  inclusion  of  f.;  without  actually  solving  for  all  the 
coefficients.  If  the  value  of* 'objective  function  is  improved 
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then  f ^  does  belong  to  the  group.  If  there  is  no  improvement 
on  the  other  hand,  the  function  f^  is  to  be  rejected.  At  the 
end  of  the  analysis  only  the  functions  that  are  of  significance 
will  remain.  At  this  stage,  we  proceed  as  in  correlation 
analysis  and  the  coefficients  can  be  evaluated  accordingly. 
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VIII.  NUMERICAL  ILLUSTRATIONS 

A.  Curve  Fitting 

In  order  to  compare  the  various  methods  we  chose 
to  use  the  data  from  "Data  Systems  and  Processing  Department 
of  Shell  Oil  Company "( 10) ,  for  which  other  curve  fitting 
techniques  have  been  applied.  The  methods  of  approximation 
are  as  follows: 

1.  power  polynomials  using  least  squares  criterion, 

2 .  self-generating  orthogonal  polynomials  using  least 
squares  criterion, 

3.  Chebyshev  polynomials  using  least  squares  criterion( 10) , 

4.  dual  linear  programming  using  Chebyshev  criterion. 

In  using  the  Chebyshev  polynomials  Tj(x)  in  the 
approximation 

i.e.  f ( x )  =  cq  +  c1T1(x)  +  c2T2(x)  +  ...  +  cnTR(x) 

the  expansion  coefficients  are  evaluated  as  follows: 


co  X  f(Xi> 

n+1 

i=0 

and 

n 

1  ^ 

c .  =  — -  y  f(x.)  t-^) 

3  n+2  U  J 

i=0 

(VIII .2 

where 

7r(i  +  1/2) 

x  - -  for  i  =  0,  1, 

n+1 

•  •  •  t  ^  ♦ 

(VIII .1) 


. 
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In  order  to  find  f(x.j),  the  given  set  of  data  has  to 

be  first  normalized  into  the  range  of  ±1.  Then  some  inter¬ 
polation  formula  is  used  to  find  f(x^) .  The  authors  of 
ref erence ( 10)  used  3rd  degree  Lagrangian  interpolation. 

If  the  variance  is  used  as  the  criterion  for  close¬ 
ness  of  fit,  all  four  methods  yield  comparable  results.  Indeed, 
since  methods  (1),  (2)  and  (3)  all  use  the  least  squares 

criterion  they  should  theoretically  have  the  same  variance. 
However,  the  use  of  Chebyshev  polynomials  might  be  expected 
to  show  a  higher  variance  due  to  the  fact  that  an  additional 
interpolation  step  is  required.  The  significant  higher 
variance  shown  by  the  power  polynomials  for  the  5th  degree 
case  is  due  to  the  calculational  difficulty  encountered  in 
solving  the  coefficient  matrix  D  of  equation  (II. 6),  which 
tends  to  be  ill-conditioned.  For  these  reasons  the  use  of 
orthogonal  polynomials  is  to  be  preferred. 

As  expected,  the  results  from  the  dual  linear 
programming  using  Chebyshev  criterion  show  a  significant 
reduction  in  the  minimal  maximum  absolute  deviation  as  com¬ 
pared  to  those  of  the  least  squares  analysis.  But  a  somewhat 
higher  variance  is  obtained. 

We  might  also  modify  the  least  squares  technique 
so  as  to  approach  the  Chebyshev  criterion,  by  using  larger 
weighting  factors  on  those  points  which  show  greater  deviations. 
When  these  weighting  factors  were  made  proportional  to  tiioir 
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TABLE  1 


Comparison  of  Results  in  Curve  Fitting 

by  Various  Methods 


Degree  of 
Approximation 

4th 

5  th 

Method  of 
Approximation 

Variance 

Max.  Abs. 
Deviation 

Variance 

Max .  Abs . 
Deviation 

Power 

Polynomial 

2.5756xl0~4 

0.0384 

1 . 3125x10 #4 

0.0240 

Self-generat¬ 
ing  Orthogonal 
Polynomials 

2  .6945x10" 4 

0.0383 

0.6831xl0"4 

0.0182 

Chebyshev 

Polynomials 

3.0936xl0~4 

0.383 

0 .87  55x10  *4 

0.0274 

Dual  Linear 
Programming 

3.9436xl0"4 

0.0233 

1 . 1190xlo"4 

0.0120 
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TABLE  2 

Effect  of  Modification  of  Weighting  Factors 

on  Maximum  Absolute  Deviations 


Degree  Maximum  Absolution  Deviations  on 

of 


Approximation 

3rd 

4th 

5  th 

6th 

(i) 

0.1041 

0.0383 

0.0274 

0.0183 

(ii) 

0.0909 

0.0383 

0.0182 

0.0110 

(iii) 

0.0579 

0.0266 

0.0160 

0.0083 

where 

(i)  from  Chebyshev  Polynomials ( 5) 

(ii)  from  self-generating  orthogonal  polynomials 

(iii)  linear  modification  of  weighting  factors  of  (ii) 

with  range  1  ^  <  5 
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absolute  deviations  some  reduction  in  the  maximum  absolute 
deviation  was  achieved  as  indicated  in  Table  2.  However, 
even  with  continuous  application  of  the  linear  modification 
procedure  these  maximum  absolute  deviations  did  not  reach  the 
true  value  obtained  by  the  dual  linear  programming. 

B .  Surface  Fitting 


In  the  case  of  two-dimensional  surface  fitting  the 


data  are  also  obtained  from  "Data  Systems  and  Processing 
Department  of  Shell  Oil  Company "( 10) ,  in  which  the  authors 
used  Chebyshev  polynomials  as  the  approximating  polynomials, 

i  .e . 


m 


n 


(VIII .3) 


i=0  j=0 


where  the  expansion  coefficients  are  evaluated  as  follows: 


n 


m 


(n+1) (m+1) 


a=0 


b=0 


but 


c 


oo 


and 


c 


with 


77 (a  4-  1/2) 


x 


a 


n+1 


Again  f(xa,  y^)  were  obtained  from  the  given  data  by  3rd  degree 
Lagrangian  interpolation . 

The  results  of  surface  fitting  with  a  4th  by  4t;h 
degree  polynomial  as  the  approximating  function  and  using  the 
least  squares  criterion,  are  3hown  in  Table  3.  As  in  the 
case  of  curve  fitting  the  self-generating  orthogonal  poly¬ 
nomials  gave  better  results  than  those  by  the  Chebyshev  poly¬ 
nomials,  This  is  probably  also  due  to  the  fact  that  the 
interpolation  of  data  is  not  necessary  for  the  former.  The  use  o 
power  polynomials  for  surface  fitting  is  not  practicable  because 
of  the  inherent  computational  difficulty. 
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TABLE  3 

Surface  Fitting  by  Least  Squares  Criterion 


n  m 

f(x,y)  =  ^  y  Pi(x)  Qj (y) 

i=0  j=0 

at  n  =  m  =  4 


Self -Generating 


X 

Y 

f  (x,  y) 
given 

Chebyshev 

Polynomial 

Orthogonal 

Polynomial 

0.2 

1.10 

0.3560 

0.3660 

0.3583 

1.2 

1.1 

3.1800 

3.2323 

3.1589 

0.4 

1.2 

0.5900 

0.585? 

0.5808 

1.4 

1.2 

2.6800 

2 .7043 

2.7043 

0.5 

1.3 

0.6050 

0.5976 

0.5952 

1.5 

1.3 

2.2300 

2.2336 

2.2316 

0.6 

1.4 

0.6000 

0.6042 

0.6041 

1.6 

1.4 

1.9100 

1.9391 

1.9296 

0.8 

1.5 

0.7100 

0.7099 

0.7048 

1.8 

1.5 

1.800 

1.7655 

1.7686 

1.0 

1.6 

0.7500 

0.7627 

0.7652 

2  .0 

1.6 

1.6300 

1.6445 

1.6300 

Mean  absolute  deviation 

0.0196 

0.0125 

Maximum 

absolute  deviation 

0.1339 

0.0538 

Maximum  percentage 

deviation 

3.95 

4.20 

C .  Evaluation  of  Benedict-Webb-- Rubin  Equ a ti on  by  Du a 1 

Linear  Programming 


If  the  approximating  functions  have  to  be  of  some 
given  form,  then  the  use  of  self -generating  or  any  other 
orthogonal  polynomials  as  approximating  functions  is  no  longer 
applicable,  because  they  contain  all  the  combinations  of 
powers  up  to  xnym.  For  instance,  if  the  equation  of  state  by 
Benedict-Webb~Rubin( 11)  is  required,  i.e. 

P(T,d)  =  RTd  +  Bq  RTd2  -  AQd2  -  cQd2/T2 

+  b  RTd3  -  ad3  +  aad6  (VIII. 5) 

2  2  I  2  \  2 

+  cd  /T  |  (1  +  yd  )  >•  exp  (-yd  ) 

then  the  coefficients  to  be  determined  are  B0,  Aq,  c  ,  b,  a, 
a,  c  and  y  with  P,  T,  d  and  R  as  the  known  quantities.  They 
can  be  obtained  by  a  least  squares  analysis  with  y  being 
assumed  arbitrarily.  Unfortunately  with  every  new  value  of 
y  tried,  the  whole  process  of  the  analysis  has  to  be  repeated. 
Thus,  it  is  very  time  consuming. 

However,  if  y  is  assumed,  then  the  coefficients  in 
equation  (VIII. 5)  will  be  in  linear  combination  and  linear 
programming  can  be  applied. 

Thus,  the  problem  reduces  to  the  following*. 


minimize  A 
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subject  to 


A  +  Eyvr^2  -  Adi2  -  co  di2/Ti2 


+  bRT.d. ^  -  ad.^  +  a,d.^ 

11  1  1  i 


+  c  e 


-Yd/ 


|  (1  +  7d±  )|di2/Ti2>  (Pi  -  RTidi) 


A  -  BRT.d.2  +  Ad.2  +  c  d.2/T.2 

0X1  1  o  l  '  1 


-  bRT.d. 3  +  ad.3  -  a.d.6 

ii  l  1  i 


-  c  e 


-Ydi' 


{  (1  +  7di2)|di2/Ti2>-(Pi  -  R^di) 


where 


for  1  <  i  <  m 


(VIII. 6) 


a  ^  =  a  a 


The  dual  problem  will  be  to  maximize 


m 


X  (pi  -  RTidi}  (ui  -  vi> 


B 


i=l 


subject  to 


m 


£  («i  +  Vi 


i=l 
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and 


m 

X  RT^2  (u±  -  v±)  =  0 

i=l 

m 

7  -  di2  (ui  -  vi)  =  0 

i=l 

m 

X  -  4 A-  («i  -  VL)  =  0 

i=l 


m 

7  RTidi3  (ui  -  vi)  =  0 

i=l 


m 


i=l 


m 


i=l 


m  2  d. 2 

V  e”Tdl  (1  +  yd±2)  (ui  “  vi)  =  0  (VIII. 7} 

.  T, 

i=l  1 

It  can  be  seen  from  the  dual  problem  of  (VIII.7) 
that  each  constraint  is  equivalent  to  an  unknown  in  the  primal. 
Hence,  a  system  of  2m  inequality  equations  of  (VIII. 6)  with 
eight  unknowns  is  reduced  to  an  equivalent  onderdetermined 


vm> 
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system  of  2m  unknowns  and  eight  equations.  Now  with  every  y 

tested,  only  the  eighth  constraint  in  equation  (VIII. 7)  is 
changed  and  the  whole  system  does  not  require  repeated  cal¬ 
culation.  Hence,  a  series  of  values  of  y  could  be  tested  un¬ 
til  the  optimum  B  is  obtained. 

Having  obtained  the  optimum  base  of  the  dual  problem, 
there  arises  the  necessity  of  finding  the  equivalent  optimum 
base  of  the  primal  which  consists  of  the  equations  represented 
by  the  basic  u.j  and  v^  .  This  system  is  then  solved  by  the 
L  -  U  matrix  inversion  which  is  further  improved  by  the  method 
suggested  by  Hotelling( 12 ) . 

The  data  for  the  approximation  are  taken  from  Douslin 
et  al(l 3)  and  the  results  are  tabulated  in  Table  4,  Both  least 
squares  analysis  and  the  dual  linear  programming  technique  were 
applied  on  the  same  set  of  data*  Fortunately,  the  least  squares 
technique  was  feasible,  since  the  matrix  resulting  from  the 
normal  equation  (II. 6)  was  well-conditioned.  However,  even 
here,  the  linear  programming  technique  took  significantly  less 
computer  time.  One  of  the  main  advantages  of  the  dual  linear 
programming  method  is  its  ability  to  handle  any  given  functions 
of  multi-independent  variables  with  data  randomly  spaced. 
Furthermore,  for  each  additional  approximating  function  intro¬ 
duced  into  the  analysis,  it  yields  a  corresponding  optimum  B, 


which  reveals 
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TABLE  4 


B- 

-W-R  Equation  by  Various  Methods 

Coefficients 

B-W-R  (11) 

Linear 

Prograiaming 

Least 

Square 

Bo 

0 .45462  5xl0_1 

0 . 409053 5 3xl0-1 

0 . 4361 52 7 8x1 0_1 

Ao 

1.79894 

1.5725351 

1.7379948 

c 

o 

0.318382xl05 

0. 4350818xl05 

0.34284477xl0S 

b 

0 . 2  52033xl0~2 

0.28653194xl0”2 

0.25722946xl0”2 

a 

0.435200X10”1 

0.61719941xl0~1 

0 . 41 38 39 82x10 _1 

aa 

0.143616xl0”4 

-4 

0.16265777x10 

-4 

0.13092259x10 

c 

0.358786xl04 

0 . 47 14486 5x1 04 

0.34 59 0712 xlO4 

y 

0.0105000 

0.0100 

0.0100 

Maximum 

Absolute 

Deviation 

4.25  atm 

0.263  atm 

0.529  atm. 

Maximum 

%  Deviation 

1.2  3 

0.354 

0.211 

Maximum  absolute  deviation  and  maximum  %  deviation  do  not 
necessarily  occur  at  the  same  point. 


r.  *  -  w 
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the  information  whether  the  approximation  is  within  the  allow¬ 
able  accuracy  or  further  terms  are  required,  without  actual 
calculation  of  all  the  expansion  coefficients. 


D.  Application  of  Dual  Linear  Programming  to  Multiple 
Regression  Analysis 

Assume  that  we  want  to  find  an  approximating  function 
f(x,y)  to  fit  a  set  of  points  (x,y)  which,  in  fact,  is  generated 
arbitrarily  from  the  following  expression 


with 


is 


F(x,y)  =  3.2  +  0.7x  -  3.0y  +  1.4xy  (VIII. 8) 

2  s?  x  <  3  and  1  <  y  <  2  . 

Now  suppose  the  approximating  function  to  be  tested 
f  ( x,  y)  =  a-L  +  a2x  +  a3xy  +  a4e-x  +  a5/y  +  a6y  (VIII.  9) 


where  are  the  expansion  coefficients  to  be  determined. 

Accordingly  the  dual  problem  is 
to  maximize 

m 

X  fUiYiMui  -  VL)  =  B 
i=l 


subject  to 


m 


i=l 


(VIII .10) 


and 

m 

^  gij-(ui  -  Vi)  =  0  for  l<j<6 

i=l 
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where  represent  the  approximating  functions  in  (VIII. 9). 

The  results  of  the  analysis  are  shown  in  Table  5,  so 
that  the  final  approximating  equation  is 

f(x,y)  =  +  a2x  +  a3xy  +  a,-/y  +  a^y 

From  the  final  optimum  base  of  the  dual  problem  the  coefficients 
can  be  obtained  from  the  corresponding  equations  in  the  primal 
problem  to  be 

f(x,y)  =  3.20  +  0.70x  +  1.40xy  +  0.884  x  10_6/y  -  3.00y 

exactly  as  in  (VIII. 8)  if  the  term  (0.884  x  10”^/y)  is  neglected 
The  analysis  was  done  with  36  data  points  in  6 
seconds  in  IBM  7040.  The  advantage  of  this  approach  over  the 
existing  methods  is  that  a  direct  test  on  a  presumed  approxi¬ 
mating  function  can  be  made,  by-passing  each  time  the  evaluation 
of  coefficients  and  the  calculation  of  variance  until  the 
final  regressional  form  is  obtained.  Furthermore,  with  the 
addition  of  a  new  constraint,  i.e.  an  introduction  of  an 
extra  approximating  function,  complete  analysis  is  not  required 
for  the  evaluation  of  the  objective  function.  Thus,  in  large 
systems  the  advantage  of  this  method  may  be  very  significant. 
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TABLE  5 

Regression  Analysis 


Constraint 

Tested 

Corresponding 

Function 

Optimum 

B 

Remark 

^io 

k 

7  .70 

gil 

1 

1.75 

accepted 

gi2 

X 

0.600 

accepted 

gx3 

xy 

0.300 

accepted 

gi4 

l/e° 

0.300 

rejected 

gi5 

i/y 

0.212 

accepted 

gi6 

y 

-7 

1  x  10 

accepted 
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CONCLUSION 


1.  If  the  criterion  of  least  squares  is  to  be  used 
and  the  approximating  function  is  to  be  a  complete  set  of  power 
polynomials,  then  the  method  of  self-generating  orthogonal 
polynomials  is  recommended.  This  method,  presented  in  the 
literature  for  functions  with  one  independent  variable,  has 
been  extended  in  this  work  to  functions  with  two  independent 


variables.  Since  data  are  usually  given  in  the  form  of  rectang¬ 
ular  grid,  i.e.  in  tabular  form,  the  developed  method  can  always 
be  applied  directly,  while  the  use  of  any  other  orthogonal 
polynomials  may  require  an  additional  procedure  involving 
data  interpolation. 

2 .  If  the  absolute  minimal  maximum  deviation  or 
Cliebyshcv  criterion  is  to  be  used,  then  the  method  of  dual 
linear  programming  is  feasible  and  computationally  efficient. 

The  method  has  been  successfully  applied  to  cases  with  one  or 
two  independent  variables.  Furthermore,  the  approximating 
function  can  be  any  linear  combination  of  specified  functions 
and  the  analysis  can  be  applied  to  any  randomly  spaced  data. 

3.  Preliminary  work  indicates  that  the  above  methods 
give  good  fit  for  the  functional  values,  but  this  does  not 
ensure  that  their  derivatives,  obtained  from  the  approximating 
functions,  will  yield  equally  reliable  results.  However,  the 
linear  programming  techniques  permit  the  inclusion  of  additional 
constraints  with  respect  to  the  derivatives.  For  instance, 


constraints  of  the  form 


n 


I 


a  .g  !  ( x .  )  ^  X 

j  j  1 


1 


can  be  readily  included  in  the  primal  problem.  Further  study 
on  this  subject  is  recommended. 


4.  Due  to  its  computational  efficiency  the  dual 


linear  programming  technique  may  offer  some  calculational 
advantages  over  the  least  squares  technique  in  conventional 
regression  analysis.  However,  only  preliminary  work  in  this 
area  has  been  done  in  this  thesis  and  hence  further  study  will 
be  required. 
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APPENDIX  A 

Derivation  and  Proof  of  the  Orthogonality  of 

Recurrence  Equation  (11.12) 


Define  the  orthogonal  polynomials  as  follows 


Po(x) 

=  1 

uQ> 

P1(x) 

=  x  Po(x)  -  Po(x) 

( iL) 

P2(x) 

=  x  P1(x)  -  a2  P-^x)  -  PQ(x) 

(12> 

•  •  • 

P^x) 

=  X  Pi.^x)  -  a.  P^U)  -  Pi_2(x) 

•  O 

(li) 

Here  the 

a.  and  B.  are  numbers  constructed  to  make 
i  l 

the  ortho- 

gonality 

relationship,  as  follows,  hold 

m 

X  W  W  =  0 

(2) 

k=l 


for  i  /  j .  Such  polynomials  are  said  to  be  orthogonal  over 
the  point  set  (x^,  x2,  •  •  *  •  xm)  *  -*-s  Prove<3  by  induction 

that  this  is  indeed  true. 

Proof : 

To  find  a-j^  there  is  the  equation 


m 


k=l 


0 


(3) 
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To  determine  a ^  and  (3 ^  two  equations  are  available 

i  .e . 

m 

X  P2(xk)  Po(xk>  =  0  (4a) 

k=l 

and 


m 

y  p2(^k)  p1(x]t)  =  o  (4b> 

k=l 

Now  P  ,  P^,  P 2  by  construction  are  mutually  orthogonal 
to  one  another,  i.e.  satisfying  (2). 

Consider  P3(x)  as  the  first  case  in  the  induction, 

P3(x)  =  x  P2(x)  -  a ^P^ (x)  -  ^P  (x)  (5) 

Multiply  (5)  by  P^(x^)  and  summing  over  k  yields 
m  m 

y  p3(^)  PjtXk)  =  X  V2(\>  pj(xkj 

k=l  k=l 

m  m 

-  a3  X  P2(*k)  Pj(V  "  Ps  X  W  P2(Xk}  (6) 

k=l  k=l 

For  j  <  3  -  2  (i.e.  j  =  0),  the  last  two  terms  in  (b)  are 
zero  by  (4) . 
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Since  x  P^(x)  is  a  polynomial  of  degree  j+1  (i.e.  1st) 
in  x,  it  can  be  expressed  as  a  linear  combination  of  polynomials 
Pq(x)  *  P1(x) ,  Thus, 


m 

Y  {  xk  Pj^xk)}  =  0  (as  j  =  0)  (7) 

k=l 


Hence  P^(x)  orthogonal  to  Pq(x)  for  any  choices  cf  and 
If  j  =  2,  take 


a 


3 


m 


k=l 


*k  P2  (xk> 


(xk) 


k=l 


(8) 


then 


m 

P3^Xk)  P2^Xk^  =  °- 

k=l 


If  j  =  1,  take 

m 

Y  XkP2^xk)  Pl^xk) 

k=l 

P3  =  _ 

m 

i  -i2'**’ 

k=l 


(9) 


(10) 


xt  >: 
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then 


X,  P3(x]^  Pl(xk) 


k-1 


0 


(11) 


Thus,  P^/  P2/  P-j  »  PQ  are  orthogonal 
Now  assuming  PQ(x) ,  ...  , 

gonal  with  ...  ,  and  # 

satisfy  the  relationship  (2),  it  is 

Pi+l(x)  . 


to  one  another . 

P^(x)  are  mutually  ortho- 
,  |3^  constructed  to 
required  to  determine 


Again  multiplying  pi+i(  x^.)  by  Pj(x^)  and  summing 
over  k,  we  get 


m  m 

I  pi+l(xk>  Pj(V  =  X  W  l  Pj(vl 


k=l 


k=l 


(12) 


m 


m 


a 


i+1 


y  w  w  -  Pi+i  x  pi-i(xk>  w 


k=l 


k=l 


For  j <  i  -  1  (i.e.  j  =  0,  1,  ...,  i-2) ,  the  last  two  terms 

in  (12)  are  zero.  Also  x-P^(x)  can  be  expressed  as  a  linear 
combination  of  PQ(x)  ,  Pj+^(x)  (which  is  still  lower  than 

P ..  (  x)  )  .  Thus, 


m 


I 


P 


0 


(13) 
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Taking  j  =  i,  and 


m 


k=l 


then 


m 

I  pi+i(xk)  pi(xk>  =  0 

k=l 

In  taking  j  =  i  -  1  and 


(14) 


(15) 


i+1 


m 


k=l 


then 


Pi+1  Pi-l^xk) 


k=l 


0 


(16) 


(17) 


Consequently  pi+1(x)  is  orthogonal  to  Pq(x) ,  .  P^x) 

if  a .  and  p.  .  are  constructed  as  (14)  and  (16)  respectively. 
i+1  l+l 

This  completes  the  proof. 
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APPENDIX  B 

Construction  of  Orthogonal  Polynomials 

of  Pi(x)-  Q-j(y) 


If  the  coefficients  a  and  6j  are  constructed 

by  equations  (III. 10)  and  the  polynomials  P^(x)  and  Qj(y) 
are  defined  by  equations  (III. 2)  and  (III. 3)  respectively#  then 
the  polynomial  products  P-(x)Q.(y)  will  satisfy  the  ortho- 
gonality  relationship  (III. 8) 

Proof : 

From  equation  (III. 8)  let 
u  v 

I* 

k=l  k=l 


i(xk>  Ps(xk> 


Qjtyj.) 


QtCyp 


(i: 


Case  1: 

t  =  j  but  s  /  i,  expanding  P^x^)  it  yields 


v 


u 


0  =  y  0t2(yr)  l  {  xicpi-i^k)'rs(xk) 


1=1 


k=l 


-  °i  Pi-l(Xk>  W  "  &i  Pi-2(Xk>'PB(Xk>} 


(2) 


It  has  been  proved  in  Appendix  A,  if  s<  i  -  2,  (2f  s  0. 


With  s 


1 


1 


and 
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u 


V1  2 
)  xk  pi-hxk> 

k=l 


U 


1  K  ) 

1-1  K 


k=l 


the  second  summation  is  zero. 
Therefore ; 


(3) 


0  =  0. 

With  s  =  i  -  2  and 


u 

7,  *k  pi-l(xk>  Pi-2(xk> 
k=l 


(4) 


the  second  summation  remains  zero,  so 


0  ~  0. 


Case  2: 


If  s  =  i,  but  t  0  j,  then  by  similar  analysis  to 


Case  1,  it  can  readily  be  shown  that  with 


' 
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X  yi 

k=l 


v 


1=1 


and 


(5) 


v 


L _ ) 

1=1 


l(Yl)Qj-2(yi) 


6  = - - - - - —  (6) 


v 


1=1 


then 

0  =  0  for  j  /  t« 

Thus,  relationship  (III .8)  is  proved. 
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APPENDIX  C 


Listings  of  Computer  Programs  and  Sample  Output 

1.  Curve  Fitting  with  Least  Squares  Criterion 

using  self-generating  orthogonal  polynomials 

2.  Curve  Fitting  with  Least  Squares  Criterion 

using  power  polynomials 

3.  Curve  Fitting  with  Chebyshev  Criterion 

using  linear  programming 

4.  Surface  Fitting  with  Least  Squares  Criterion 

using  self-generating  orthogonal  polynomials 

5.  Evaluation  of  Benedict-Webb-Rubin  Equation  of  State 
by  Linear  Programming 

6.  Evaluation  of  Benedict-Webb-Rubin  Equation  of  State 
by  Least  Squares  Analysis 

7.  Multiple  Regression  Analysis  by  Linear  Programming 


^-o  i 


%  JOR 
%  i  r Mr 

$  r  H  JOH  OR  I  HOG 

%  i  « I-  i  c  r  u  p  v  t- 


lSuuui  CURVp  Pill  I i\iG  p  t  iii-  Mi- ■  km  i  I  iM  ORTmuGoNAl  kuL  i  mOHIAL 
* you 
GO 

N  O  O  F  (  K  *  MH|  I  s  | 


( 

r 

r 

C 

r 

C 

( 

c 

c 

r 

r 

r 

r 

r 

C 

r 

c 


i  u  o 

i 


2 

1  01 
•a 


Vl  PITTING  r  SF  LF-GP  NFR AT  I N6  .  f  1  jGi  >NA l.  POL»  IM  )'  I  A l  • 
vi  -  f\io .  OP  OAT  A  MO  f  's|  1  *  aMu  1  00  • 

N  =  O  P  G  K  P  P  OP  APRROX  i  -I  A  i'  I  ON  *  ■  '■  a  |  uM  in, 

a  a  i  i  /  *  -ha  i  i  i  ami-  THF  COFFp  I C I FNTS  OF  THE  ORTHOGONAL  POL 
ill)  A  k  f  Tpf  P  r  a  m  A  N  -  i  O  N  COtFFICIENTS  0  ■  *  i  0  iM  * 

WITH  A  X  (  1  )  A  IM  U  pi  A  (  J  )  • 

All)  A  R  P  i H  r  l  N  0  F  H  E  N  U  t  M  T  V  A  R 1  A  6  L  P  S • 

V  (  1  )  Amp  i  pip  DFPFiMDF.M  i  VARIABLES. 

A  SURPROgRAM  IN  C'ONVFRTlNG  Thi-  ORTHOGONAL  POL  r  in' »"!  I  AL  s  INTO 
POWFR  POI  YNOMTAI  S  |S  A I  SO  IMrLUOpO. 
r  H  A  NGF  TNPUl  AND  OUTPUT  PORMAT  |p  NPCpGGARr. 

I last  oa i a  roM T mol  rARn 

+  V1-  FOR  C  PI  A  N  G  P  OF  D  P  G  R  P  P  OF  APpROXIMAiI(  vi  ONLY# 

ZPRO  POP  NPW  S P i  OF  DATA* 

—  VP  To  PNO  I  HP  (  QMPU  F  A  i  1  OIM  • 


O I M  FNG I  ON  A (100)*  r ( 1 0  0 )  *  Axil  5)*  PA (15)*  C(15)*  pi  15*100) 
RF  Al )  (5*1)  M 

FORMA  I  (  2  X  *  1  5  ) 

FM  =  M 

RF  AO  (  5  ,  ?  )  (  a  (  1  )  *  Y  l  i  )  *  I  =  1  *M  ) 

FORMAT  (?X*  PPM*  OX*  PH.4) 

R  F  A  I )  M  *  1  )  N 
W  P  T  T  F  (  A  *  4  )  N 

FORMAT  ( l H—  *  5X*  2  6  H  f ) F G R F F  OF  APPROX]  'AT  ION  =  *  i  *, 

MA  =N+1 


F  N  =  N  A 

(  G  E N  fc  R  A  T 1 0 N  O F  ORTH 0 G O  N A  L  POL V N O M I  A L  S • 

DO  10  K  = 1  *  M 
lu  P (  1  ,<  )  =1 . 00 
1  =  2 

SI  =  o , O 0 
F  ?  -  ! )  .  o  n 
no  1  ?  K  =  1  *  M 

s  1  =  s  1  +  X  (  <  ) 


5  2  =  S  2 + 1  . 

A  A  (  I  ) =S1 /S2 
R A  (  1  )  =0 . UU 
no  l 4  K  = i *  * 


1  3 

P  (  I  *  K  ) 

=  X ( K) - 

A  A  i  1  ) 

w  p  t  r  f 

(6*14) 

1  4 

FORMA l 

(  1  HU  * 

5  X  * 

39HC0FP P IC I EN i S  Or  QRT HOGONAL 

WRITE 

(  A  *  1  5  ) 

1  5 

FORMAT 

(  1  H  * 

HA  ♦ 

H  H  A  a  (  I  )  *  MX*  5  H  P  X  (  I  )  ) 

WR  I  IP 

( P  *  ip) 

AX  (  1 

)  *  P  X  (  1  ) 

1  A 

FORMA  1 

(  1  H  * 

?  MX 

*  PI s.H)  ) 

OO  2  U 

I  =  H  *  M  A 

T  A  =  T  -  1 


i  I  AL  • 


. 

•  t 


• 

(  . 

•  •  * 

• 

•  « 

(  .  )  - 

. 

• 

. 

• 

(  • 

•  X 

» 


r  . 

. 

/L~bc 


M  =  n .  u> 
c2  =  n.  ) 
f  3  =  o .  o  o 

S4  =  0.  )0 

r>o  IQ  K  =  1  , M 

Sl=b1+X(K)*P(  f-l  »  K  )  *•*  2 
b2=b2+P  (  T  —  1  %  K  )  **2 
S  X  =  b 3  +  X  (  <  )  *P  (  I  -  i  »  K  )  *P  (  I -2  9  < ) 

19  b4=b4+P ( I— 29<)**2 

A  X  (  I  A  )  =  b  1  /  b  2 
R  X  (  T  A  )  =  b  3  /  b  4 

W  P  TIP  1X921)  A  X  l  J  A  )  ,  hxi  |  A) 

2  1  F  O  P  M  A  I  (IN  ,  2  (  4  X  ,  F  •  .  H  )  ) 

no  2  U  K  =  |  ,  M 

P(  I  *K  )  =  (  X  (  K  I  —  A  X  (  I  A  )  )  *  P  (  1-1  9  K  ;  -  R  X  (  |  A  Ml-  P  (  |  —  2  »  •<  ) 

T  r  A  L  CUL  A ‘I  I  OIM  OF  FXPA.MbfOIN  (  OF  F  F  J T  1  FlMTb. 

W  P  T  j  F  (  A  *  2  2  ) 

FORMA  I  (1HU*  '-'X*  I  H  J  9  H  A  9  4HC  (  I  )  ) 

1)0  2  3  J=  ]  9  M  A 
( I =0 . 00 

r?=o . uu 

DO  2  4  k  = i *  M 

n.  =  ci  +  y(K)*p(j,< ) 

24  r 2 =02+P ( J • K ) **2 

r  (  j)=n  /c  2 

WR  I  TF  (  X.  ,  2M  J  9  r  (  J  ) 

2  u  FORMA  I  (  i  H  «  4X»  12*  X  9  F  ]  *  H  ) 

r  ffmok'  of  approx  i  mai  ion. 

WPT  ib  (  X ,  2  6  ) 

2X  FORMA  I  (  I  FlU  9  Xa»  ]  FI  X  9  8X9  I  H  1  9  X  X  9  X  HC  A  L  #  1*  3  X  9  3  H  *  » F  v/  ) 

no  30  I =1 , M 

W  =  0  .  U  I 

no  A  I  J  =  i  »  IN  A 
A  1  W=  W  +  C  ( J ) *p ( J  »  I  ) 
l)=  ARb  (  tf-Y  (  1  )  ) 

"AO  W P  T  I  F  (  X  *  b  2  )  X  (  )  )  *  2  (  T  )  »  W  9  F) 

3  2  F  ( )  R  M  A  1  (  1  FI  9  3X»  F  X  .  2  9  ^  X  9  F 6>  •  2  9  3  X  •  F  *  «  4  9  4X9  F  8  •  4  ) 

r  rONVFRS  TON  OF  ORTHOGONAI  ro  POWFR  POl  YNOMIALS* 

PI  1  90=1.00 
P  (  2 , 2  )  =  1  .  0  0 
P  (  2  9  1  )  =  —  A  X  l  1  ) 

P  (  3  9  3 ) = 1  , 00 

P(  3  9  2 ) =  P (  2 ♦ 1  )-Aa(2)*P(2*2) 

P(  |  )  =-p  (  2  9  I  )*AX(2)-fX(2) 

DO  4  0  I  =4  9  IN  A 
I  A  =  I  —  1 
JA=  |-2 

p  (  j  ,  r  )  =  p  (  1  a  .  1  a  > 

P  (  T  ,  T  A  )  =  P  (  1  A  9  J  A  )  -  A  X  (  I  A  )  *  P  (  1  A  9  i  A  ) 

UO  4  1  J  =  2  9  J A 
JH=  |  -J 

41  P(  |  9  J  R  )  —  P  l  i  A  9  Jp  — 1  )-AX(  1APKI  i  A  9  J  fO  —  P  A  (  1A)*P(JA9Jp) 

4  0  P  (  1  9  1  )  =  -  A  X  l  J  A  )  *  P  l  I  A  9  I  )  -  H  X  l  I  A  )  *  P  (  J  A  9  1  ) 


.  . 


*  (  ^  *  •  - 

• 

1  • 

•  ( 

.  u 

• 

• 

•  *  ) 

•  * 

.  • 

. 

• 

r  . 

•  i 

.  I 

•  i 

9 

>  <  M  1 

• 

. 

• 

(  •  ’  )  U  'M 

•  •  I  •  (  > 

.  ,  *  .  .  . 

.'tfr.riu 
.fsfC.C  ) 

-s (  r . C  »  q 
.  r  r  r.i  )  q 

•  I  (  )  '  t  V.  -  {  »  ^  )  -t 

(  .  )  -  -  1  •  I 

#  ^  1 

(  •  I  •  t  )  »-  (  T  .  r  >  4 

.  ft.  14:1  *  .  T  I  -I 

.  V 

=  '  r  *  I  )  *> 


2.-6  3 


WR  T  TP 

(6*42) 

42 

forma i 

(  1  H-  , 

2  X  *  4  4  H  F  T  N  A  L  f 

WR  T  1  F 

(6*43) 

4  3 

FORMA  1 

l  i  H  o  * 

PX*  1 H ! ,  lux* 

IX)  6  u 

1  =  1  *  IMA 

AX (  I  )  = 

u  .00 

no  7  1 

J=  I  *  IMA 

2  I 

AX ( J ) = 

A  X  l  I  )  +( 

(  J  )  *  R  (  J  *  1  ) 

6  0 

WR  I  (  h 

(6*62) 

1  *  AX  (  1  ) 

6  2 

FORMA i 

(  1  H  * 

/X*  12*  hx*  hi 

C 

C  H  F  (  K 

Oh  T  Mh 

WR  I  1  F 

(6*26) 

W  =  .)  .  a) 

no  60 

T  =  1  *  M 

/  =  ) . 

IF  (  X  ( 

T  ) )  61  * 

6  2*61. 

62 

6  =  AX (  1 

) 

GO  TO 

6  3 

6  1 

1)0  6  4 

<1 

7 

* 

— 1 

II 

*> 

JA= J-1 

64 

2  =  2  + AX ( J ) *X (  | 

)  **JA 

63 

D  =  A  R  s  (  /  —  Y  (  1  )  ) 

w = w+n*n 

6  U 

w  R  T  I  F 

(6*32) 

X  (  1  )  *  Y  (  t  )  ,  /  , 

V AR=W/ ( FM-FN ) 

W  P  1  IF 

(6,66) 

V  A  R 

66 

format 

(  1  H-  * 

2  X  *  ]1  H  V  A  R  i  A  N  C 

P  1-  AO  ( 

^  *66  )  I 

D 

66 

FORMA  1 

( 2  X  *  12) 

h  '  w  F  R  R  o  j_  V  f\j  c  i  ■  /  |  A  L  2 


4HA ( i )  ) 


ft  ) 


I) 


. 


IF  (TP)  6 
6  7  COIM  I  I  i\!Uh 
t  N I) 


*  I  OO  » L  0] 


CUR  Vt 


1  .  oouo 

1  .414? 

1  .7321 

2  . 

2,23  • 
2.44 
2.64  5 ft 
2.8284 

3.1623 

3.3166 

3,4641 

. 

.  • 

3  .  ft  7  3  0 

4  . 

4.1231 
4 . 2  4  2  6 
4.  4  6  6< 


0 


I  •  '  ) 

.  V  .  . 

<  ■  0  ■  ) 

0  0 

0 

.  f  a- 1 

.  .  • 

. 


»' ,  -  r 

•  r 


•  A  ) 

• 

f  >  *  •  i  ) 

4“64 


2  1  . 

7  ]  . 

9  ?  .  )  0 

2  4  •  U  U 
•  uO 

4 

1 

5 

-1 


4.4721 

4  .  *  8  2  ft 
4  •  ft 

4. 7958 

. 

. 


(SAMPLE  Ql  JTPU  I  ) 


I ) E G R E 1  Or-  APPROXIMATION  =  5 


COEFFICIENTS  Or 
AX  (  I  ) 

•  i  30  ■  . 

0.14  F  +  2 


ORTHOGONAL  POL> 
B  X  (  I  ) 

0. 

.  f  +  ; 


.  l  J  +02  0.414  F  +  0 i 

.  i  3  P+0  2  0 . 396  ) F  + 

U  .  1  2  9  •  .  ) . 


I 

C  (  1  ) 

1 

0  • 

442  9  35  58  F  +  0 ] 

2 

0. 

1  9  4  /  3  7  7  7  h 

-00 

4 

0  • 

592900865 

-02 

4 

u. 

2 1096 4  2  U  t 

-0  3 

6 

• 

|  44  2  6  04  7  F 

-  0  4 

6 

• 

]  1  1  0  1  6  4  0  F 

-0  5 

X 

Y 

r  AL  .  Y 

D 

1.00 

1.00 

0 

'4  .  )0 

1.41 

1 . 4959 

0 

4  ,  UU 

1.74 

1  .  7  2  0  0 

0 

4.  ,)0 

2.00 

1 .9991 

0 

-s 

• 

c 

2.24 

2 . 2  4  2  8 

0 

6  .  JU 

2.45 

2*458 9 

0 

7. 

2.69 

2.6939 

« .  uO 

2.84 

2.8330 

0 

9.  0 

4 . 00 

3  •  )  )  0  4 

0 

1  ) .  ) 

4.16 

-i  .  i  5  8  9 

0 

11.(0 

4.42 

•<.4104 

12.  10 

5.45 

3*45  ;  4 

0 

14.00 

5.61 

5.6  50 

1  4  .  .10 

.  7  4 

5 .  5  8  3 

. 

4.87 

i  •  8  7  2  5 

0 

16.  >0 

4.00 

4  •  0 ' 1 2  4 

0 

.  I A  L  5  • 


tV 

.01 

. 

. 

.0009 

.0067 

•  0094 

•  0046 

. 

•  0  0  14 

.  1 

.  006  5 

.  0  0  5  6 

. 
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•  . 

. 

• 

.  r 

.  r 

. 

. 

-  • 

• 

-  •  ; 
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. 

17,  U 

4,1  7 

4  •  !  7  /  9 

0 . 00 4 H 

. 

A.  7 4 

4  •  2  4  8  '3 

0  .  5  0^7 

4  •  3  6 

4  •  3  6  >  7 

. 

2  0 . 0  U 

4  •  4  "7 

4  •  4  ~  4  •• 

- .  00 2  7 

2  1  •  u 

4  •  5 « 

4  •  8  i  : 

j  • 

2  /  .  j  0 

4.69 

4*6858 

• 

4  •  3  0 

4.7894 

0*0 0 6  4 

74. 

4  .  9  : ) 

4  •  «  P  5  0 

i  . 0040 

?^,on 

6  * 

8  .  106  7 

0 .0067 

F  T  N  A  L 

COF  FF  |f I F N 1 S  OF  i HF  POWF p  POLYNOM I 

T 

A  (  I  ) 

1 

0.  6692618E+00 

7 

0 . 4862464,/  E-00 

•a 

-0. 41823675  F - 0 1 

4 

•  7646  7  i  6 6 F  —  0  2 

6 

-0 •  8  6  5  8  6  7  0  8  E  —  0  4 

6 

0.11101 6  4  0 8—85 

X 

Y 

r  a  l  .  v 

I) F  \7 

1.00 

1  .  0 

1 .0169 

0.0169 

2.  UO 

1.41 

1  .  3  P  6  P 

0.0189 

7.00 

1.73 

1.7700 

0.0121 

4.  JO 

•  o  u 

1  .  9991 

0 .  u  u  0  9 

5  •  0  0 

2.24 

2.2428 

0. 0067 

6  • 

2.45 

2.4689 

. 

7.65 

7.6689 

0 . 

«  • 

7.83 

7.98  -s  () 

0. 0046 

7,70 

3  .  ■  0 

4  .  -404 

0.0003 

1  C  .  )0 

3.16 

0.1689 

9 .0004 

11.00 

8.3  7 

3.31 o 9 

0.0068 

1  7 . 0  0 

3.4  6 

*  •  4  ’ 

0 .  006  6 

1  8  .  J  U 

3.61 
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. 

14.00 

3.74 

3  .  7  1  8  3 

. 

1  8  ,  O  0 

8.87 

3.9  776 

0 . 0006 

16.00 

4.  10 

4.0026 

0.0 0 7  8 

17.00 

4.12 

4.1279 

0. 0048 

18.00 

4.24 

4.2483 

. 

■. 

4.36 

4  .  '  6  3  9 

. 

?  0 .00 

4.4/ 

4.4748 

. 

71.  0 

4.88 

4.681 7 

0. 0009 

7  7.9  0 

4  •  6  P 

4 . 68  6  p 

0 . 0046 

7-2  .  0 

4.80 

4.7894 

0 . 0  0  4 

7  4. 0  0 

4.  PO 

4.8980 

0. 0040 

7  8  .  0 

5.00 
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( )  .  '06  ' 

VAP  I  ANCF 
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LEAST  SoUapf  ANALYSIS  USING  POWFP  POLYNOMIALS  A 6 
APPROX IMAT  I NG  FUNCT  I  ONS  . 

H  ANALYSIS  T  S  f )  0 N E  1 Y  T  H E  L - U  M A T R I X  I N  • 

i  Hi  EXPANSION  COE  EE  I  C  I  E"  NTS  IN  ASCENOiNG  POWER  OF  X 
ARE  GIVEN  UNGER  E  INAL  SOLUTIONS. 


(SAMPLE  OUTPUT) 


i  r  •  I  FAST  SOUAPc  ANA  |  YSTS  -  ■  PO|  i  NOM I  A|  S. 


OFGRPF  of  approximation  =  5 

FINAL  SOLUTIONS. 

0 . 56745699F+00 
0 .4R7P^1 SI F  — 0  0 
-0 .41 742A23F-01 
0 . 763976 26F-0? 

-('.PAT  7;fi  77F-04 
0.11  )31  734F-01 


X 

Y 

1  •  )  0 1 

U  .  1  0  0  (J  0  i )  0  F  4  0  1 

2  •  0  9  0  u 

U  ,  L  4  1  4  2  0  U  U  E  +  0  1. 

0,  1.732  1  OUOF+O  1 

4.  )  0  0 

0 . 2  0  0  U  0 0  0 0  F  +  0 1 

5  •  u  u  0  u 

.  .  ?  3  1 

A, OOOO 

0.244O801  !  f"  +  1 

/  .  0  1 0 

.26458  4-  1 

8,000 0 

) ,  2  8  7  8  4  '  F  4- 

Q  .  0  ’  n  ' 

.3  F+ 0 1 

. 

0 , 7  1  A  9  3  0  -  0  F  +  1 

1  1  .  ■  ■ 

0.331  A  A  0  )  F  4-  0  1 

1  2.0 0 0  0 

0 , 3  4  A  4  T  0  0  0  F  4-01 

i  3 . 

0,3d 

1  4  •  u  u  0  0 

0,37417 OOOF+O i 

1  A , 0 000 

0.38  7  300  + 

1  A  .  . 

0,400000  !  1  F 

1  7.0  0  00 

0 . 41  2  3 1  0  0  F  4-  0  1 

18.00 

U  ,  4  7  47  AO  E!  OF  4-0  1 

1 9,0000 

."1.436890-  F  4-  1 

2  U , 0 UO 

0. 44  7  2  1  0  0  IF 4-  1 

2  1  • 

0.46*  6  4-01 

92.000) 

0.4 A 0040  'F+01 

9  3. 

0.47Q6Q000F+O1 

9  4  ,  U  J  0  U 

0 . 4  P  Q 9 0 0  OOF  +  01 

. 

0 . 50000000F+ 

C  A I  .  Y 

0  .  10] 60999F  +  0] 
0 • 1 3959 1 9 7  F  +  0 1 
0 • ] 7198154F+01 
0.19  8  2 4F4 

.  •  4  E+0 

,26  + 

0  . ?  P  3  30  6  5 I F4- 0 ! 

0.316Q363<  ■  - 

,111  i  a n 4 1  +ni 

.  ,i  -  ...  + 

0.3601 61 91 F  +  ) 

.  -  4 

...  302  f+ 

0.4  63492F4 

0.41 327 0 1 1F+0 1 
0.47543  '  F+ 
.4*3  1  84 

.  4  + 

.  4  A  <•>  V  A^4f  +  0  1 

.  4  fl  I  4  A  F+ 

.4  0 4  •  ^  F  + 

0.60239PH7E+  )  I 


DEV. 

0 . 1  6<  -  1 

0.182  •  - ■  1 

0.1  ?  2  8  4  5  7  7  F  -  0 1 
0.1  l 1 7  6  0  2  0  F —02 
0. 644493 1 0( —  2 
) . 9 1 8 8  •  01 E-02 
•  99  31  1 —  1  .- 

9.466  8 0  7 67F—  v 

r\ '  66  7  P  a  5  7 1 - n  3 

n  .  7  q  46  0  7  p  p  c  _  r\  9 

n  .  4  Q  9  6  4  o  ^  i  c - n  9 

■  I  .  K  4  7  6  3  P  4  F  —  '  / 
)  *  3  -  2 

0 .11531  -02 

0.2530  I P 7  4  F - 0  2 
0.6349206 0 f - o 2 
.960111 6 • ■ “02 

. 

0,1/  4444  “01 

.  -  6  - 

.  *4  -  ■  • 

m  7  7  6  3  4  3  P  6F- 
.  7  9143  6  7 4F-0 ? 
.  | 4P  -01 

0 . 2  3  9  P  H  7  2  4*  -  1 


VAR  1  ANTE 


0.13 L25054E-0 4 


<« . 

2-4  t 


rUoVf  FTTTTNG  pY  CHFBYSHFV  CR  I  TFR  ION  USING  POWFR  POLYNOMIALS* 


PART  I 
FOR  FATH 


-  LINEAR  PROGRAMMING  BY  rI'PLEv 
POWER  TERM  I NT  ROOUCED  AN  OPTIMAL 
P  0  N  0  I  ,  .  • 


Tl  I  , . 

11 A L U E  AND  m 


( SAMPLE  Ou i PUT  ) 

CHEBYSHEV  APPROXI  A I  ION  R’t  DUAL  LINEAR  PROGRAM  I  -  • 


J  h c  T  |VT 

VA  L  UE 

=  • 

2u0u0G0GF+01 

1  NOt P • 

Vector 

2  5 

0  •  5  Of 

f  N 1 )  F  P  . 

VEC TOP 

2  6 

•  50 

OpT  T  MA| 

\/A  (  UF 

0 . 

•r^-^UE-00 

T  NOFP . 

VFC  TOP 

Q 

3 , 5C  F-C 

TNOFP, 

\/rr  TOP 

?  6 

•  4334  I  -  T 

I  NHFP, 

VICTOR 

ft  0 

.  6  7 

ft  ft  * 

ft  ft 

OPT  I  M A  L 

VAL  Ue 

• 

1  4  4  7  G  8  4 ~  —  0  0 

I  N  l )  E  P  . 

V  F  C  TOR 

4 

G.41558442Ft00 

ID  f  P. 

VFCTOR 

26 

0.295454! 4 

INDEP. 

VFCTOR 

42 

• .  4  •  46  t( 

I  NO tP. 

VFCTOR 

25 

Q ,8441 5  5  F  6f —  ^ 

ft  ft  ft 

ft  ft 

OPT  IMA! 

VALUE- 

ii 

o 

• 

4g HQ862 4 F  —  01 

,  DFP. 

VECTOR 

•37758444  - 

I  N  D  F  P  . 

VP (  TQR 

2  6 

.24 

T  N  0  r-  P  . 

V  pp  TOP 

4  6 

. 

INDHP. 

V  E  (  TOP 

2  0 

>•122416  - 1 

I  MO f  P. 

VFC TOP 

50 

0 .4492  49  5 1  P-0  1 

$  ft  ft 

ft  ft 

0  P  i  I  M  A  L 

V  A  LUE 

=  u  • 

2326Ibl6F-01 

1  lM  1")  E  P  • 

Ve(  Tup 

2 

0 *35876345 F - 0  0 

i NUtP. 

VFC TOP 

26 

•  2  -  ■ 

,  HEP. 

VFCTOR 

42 

. 

I  N  D  F  P  . 

V  F  r TOP 

1  5 

.  4948F- 

T  N  D  F  P  . 

V PC  TOP 

4  7 

.  -  ' 

I  N  D I  P  . 

VPC  Top 

2  F 

ft  ft  ft 

ft  ft 

• 

•  M  • 

. 

.  U  •  • 

. 

2-6 


The 


HART  T I  -  SOLUTION  OF  PRIMAL  BASIS  PY  L-U 
£  a H  A N S  i  0 N  COt-  FF  I  C  I  t.N  i  8  A  R  r-  (j  I  V  r  N  U N D r  R  i  H F 
wjTh  i  HF  FIRST  VALUE-  A r  THE  ABSOLUTE  MINIMAL 
1  )F  V  1  A  r  I  ON  XI  T  ]  .  • 


’  *  A  T  R  I  X  I  N  V  F  R  “ 

—  t  M  a  I  C  n,  I  »  '  "r  T  ~ 

r  i  J v» A L  cuLui  [l 

*  »  A  v/  T  '  *  '  !  •  * 

AX  I  ,'j 


i 


1 AUOOl  COFFF  Id  FNTS  FOR  POWF  R  POL  Y  \  )M  I  A L  FROM  LIU  A R  PR(  .  . 

FINAL  SOLUTIONS. 

0 . 2  x  ?  6  1  6  ?  2  F  - 1 M 
0.6U3P1616F+00 

o  .  4476  zqjt+C-O  ; 


-  . 

. 

-0 . 1 93021 

27P-04 

X 

F  (  X  ) 

O  A  L  .  F 

DFV 

c 

c 

• 

1 — 

1  .  oo" 

1.0233 

•  2  : 

a  i  )  v.; 

1  .4142 

1  .  8909 

>  •  ; 

• 

o 

o 

1  .  • 

t  .  71 

.  0 

4  .  U0 

• 

1  .997] 

0 . 0  0  7  9 

5.00 

2.2361 

2.247 o 

0 .01  14 

6 . 00 

.  44 

2.4699 

7.00 

2.64^8 

2 • 66 71 

0.0233 

8  •  UO 

2.8284 

2 • 84°4 

• 

9.00 

3  • 

'.0147 

•  147 

1  0 . 0  0 

• 

3  •  1  6  8  5 

• 

1  1  .00 

■  316  6 

.  3  1  3  6 

. 

1  -  .  0  o 

8 . 464 1 

3 . 4  :  4 

0.011“ 

18.0 

•  . 

3  .  T  8  7  1 

. 

1  4.00 

3.741  7 

.  .  •  2 

.  G , 

1  5 . 0  J 

.84  9  ’ 

.0; 

16.00 

4  •  0  0  1 0 

3 . 9" 3  8 

. 

1  7  •  0  0 

4.17 

4.1081 

0.0:50 

18.00 

4.1  - 

4  •  2  3  r  9 

. 

19.00 

4  .  3  8p' 

4 .  3  61  - 

. 

2^.00 

4.4  r, 

4.4846 

• 

7  1  .  0 

4  • 

4.6026 

.  . 

2  2.0  0 

4.6904 

4 . 7  1  3  7 

. 

2  8  .  0  0 

4.79^8 

4 . 8  1  “2 

. 

7  4.00 

4  #  f 

4.9041 

. 

2  a  .  o  1 

8  .  0  0  ■' 

4.9-47 

. 

V  A  p  T  A  NTH  =  0. 

3049 88  F  ^ F - 0  5 

*  *  *  * 


• 

• 1 

.  ‘ 

1 

•  - 

• 

. 

. 

• 

• 

. 

. 

. 

$  JOR 
$ T 

J  B  JOB 
$  f  R  F  T  C 
C 
C 
C. 
r 
r 
c 

c 

c. 


1  U'J 

1 

? 


4 

4 

C 

1  0 
c 
c 


1  1 


)  4 


I  2 


I  4 


1  a 


]  5  0  ^  0  3  o  U  R  F  A  C  E  F  I  T  T  I  N  G 

2*250 

LEUNG  GO 

Surface  node ck  * nol i ' t 


£.  ~  O  - 


ri 0 G 0 N A L  t'CL  Y;\tC‘  •'  I  A L  S 


SURFA  riii  ING  WITH  SFLF—  ATI  nIG  orthogonal  poly  al  • 

IN  POl  I  -  I  . 

N  = 0  A.  KU=NO.  OF  POINTS  ALONG  THF  X-AXIS. 

Mel  iRPF  IN  Y.  LV=NO.  Y-AXIS. 

IF  T MF  |  AST  DATA  rONTROl  FARO  I S  POSlTjVF*  THF  PROGRAM  P 
A  N p W  S F  T  r 1  i  > i  T  a  O  F n 11  R  vv  I  F  IF  4  •  I  N A T F  '  • 

I  -  I)  THF  FORMATS  W(  I 


l)I  M  F  NS  [  ON  a  (  2  0  )  *  Y  i  2  ^  )  *  2.(20*20)  *  GY  (  ^  0  ) 
DOUBLE  P  R  E  (  I S 1 0 N  P  ( 1 0  *  2  0  )  *  0(10*20)*  A X ( 1 0 ) 
P.  Y  (  1  0  )  *  A  (  10*1  0  )  ♦  C  (  1  0  *  1  0  ) 

READ  (5,1  )  N  *  VI  *  <  U  *  L  V 
forma  r  ( 4 ( 2X  *  12)) 

w  R  T  T  F  (  A  »  2  )  N  ,  M 

FORMAT  ( 1 H—  *  5X>  2AHOFGPFF  Oh  APPROXIMATION 

P  F  A  0  (5,2)  ( X ( < )  *  K  = 1  , K U ) 

R  F  A  D  (5,4)  (0/(1),  |  =  l  ,  L  V  ) 

FORMAT  ( 2  X  *  F  5 . 3  ) 

R  F  A  n  (5,4)  (  ( Z  l  <  >  L  )  *  <  *  1  *  K  U )  *  |  = 1  » I  V ) 

FORMA  1  ( 1  1  FA . 2  ) 

GENERATION  OF  ORTHOGONAL  -  Li  IALS. 

DO  10  <= 1  *  KU 
P  (  i  ,  «  )  =  1 . 0  u 

IN  CA5F  OF  Thf  INVFRSF  OF  Y  IS  RFQUJRED  IN 
L  F  T  Y(L)=1.00/GY(L) 

no  11  !  =1  *  LV 
Y ( L ) =GY ( L ) 

0(1,1  )  =  1  .  )  0 
l  =  2 
J=  2 

SI =0.00 
S  2  =  0  •  0  0 
DU  12  K  =  i *  Ku 
S  1  =  S  3  +  X  (  K  ) 

.5  2  =  52+  I  .0 
A  X  (  1  )  =  S  i  /  S  2 
B  X (  1  )  =  .00 
DO  12  K  =  1  ♦ KU 
P (  I  , K ) = A ( < ) —AX ( 1  ) 

1  1=0.  )U 
12  =  0.  '  O 

no  14  (  = 1 , LV 

T  1  =  T  1  +  Y  (  L  ) 

r 2  =  f 2  + I  .00 

A  Y  l  1  )  =  T  1  /  T  2 

k  Y  (  I  )  =  u  .  u  0 

1)0  15  L  =  i  *  LV 

0  (  J  *L  )  =  Y  (  L.  )  -  A  Y  (  I  ) 

WM  [  lb  (  A  *  I  b  ) 


F  X ( 1 0 )  *  A  Y (10)* 


-  ,  12*  2  h  X  *  T  2 ) 


HF  APPROXIMATION 


.  t 

• 

f 

r  # 

•  * 

(  i » 

< 's  •  A  J 

0 

)  (v  . 

f  X  ) 

\  # 

•  *  =  • 

•  I  -  (  J  »  I 

I  *  -  I 

(  i»  r  i 

*  I  - 

.  i  = 

)  • ;  r  (  y  •  r 

.  +Os' 

'  •  I 

(  I#  K 

(  •  • 


c.  /  u 


1  6  FOR  M  A  i  (  1  H  0  *  5  X  *  39HCOt)*FICI  .  GO  N  A  L  POLYNOMIAL^#  ) 

w  R  I  l  F  (6*17) 

17  FORMA i  ( 1  H  * 8 a  *  c H  A  X (  |  )  »  j  3  a  *  5FiRX(I)*I3X*  HA'(J)*13a*  " H  r  { ( J )  ) 
wR  [If  (6*18)  A  X ( 1  )  »  B  X ( ]  ) *  A  Y ( 1  ) ,  -  7(1) 

18  FORMAT  (  1  H  *  4  (  -<  X  *  F  J  *“>  •  8  )  ) 

NA=N+1 

M  A  =  M  +  1 

no  70  T  =  3  *  N A 

i  A  =  [  -  1 

SI = 0 • o  0 
S  2  =  o  •  j  o 
S 1=0*00 
S  4  =  0 . 0  0 

no  19  K  =  1  *  KU 

S 1 =  S ] +X ( K ) *R (  1-1  *  IF  )  *  *  2 

S  2  =  S  7  +  R (  1-1  *K ) **2 

S3=S3  +  X(K)*P(  1-1  *  K ) *  P (  I - 2  *  K ) 

19  S4  =  64  +  P ( I -2  * K ) **2 
AX(IA)=S1/S2 

RX (  T  A ) =  S 3/S4 

WR I T  F  (6*21)  A X (  I  A )  *  R  X  (  I  A  ) 

71  FORMAT  ( 1 H  *  7 ( 3 X ♦  F 1 F  #  8 )  ) 

DO  2  0  K  =  1  *  KU 

20  R(  i »K)  =  (a(K)-AX(  i A )  ) *  R (  I - 1  *  K ) - P  X { I A ) *  h (  I - 2  *  K ) 

1)0  2  2  J=3»'viA 

J  A  =  J  —  1 
T  1  =  0 . 0  0 
1 2=0*0  ' 

1  3  =  0  .  )0 

T  4  =  0  .  J  ) 

no  7  3  I  =1  ♦  L  v 

T 1 =  T 1 + Y ( L ) *0 ( J— 1  *  L ) **7 

T  2  =  T2-FO  (  J-l  *  L ) **2 

T  3  =  T 3+Y ( L ) *0 ( J-l *  L ) *0 ( J-2  *  L ) 

2  3  T4  =  T4  +  0  (  J-2  *  L  )  ■**2 

AT  (  J  A  )  =  T  i  /  T  2 
RT  i  J  A ) =  T 3  /  T  4 

w  R i 7 1  (6*24)  A 1  ( J  A ) ♦  R  Y ( J  A ) 

24  FORMA'*  ( 1 h  »  3  6  X  *  2  (  3  x  *  *6)) 

DO  7  2  L  =  1  *  LV 

22  0(J*L  )  =  (  Y  (  L  )  —  A  Y  (  J  A  )  )  *Q  (  J-l  *  L  >  —  ^  t  (ja)*-,i(J-/*L) 

C  A  LCul.  A  T  1  ON  OF  F XP  AN S  1  ON  ("OFF  F  1 1  1 1  NT  . 

WRITF  (6*250) 

2  r  u  FORMAT  ilHJ*  '^H  J  PANS  ION  OOFFF  IT  IFN  .  ) 

DO  2  6  I  =  1  *  NA 
DO  25  J=1 *MA 
61=0.0., 

S  2  =  0 . 0  0 
DO  26  IF  =  1  *  Ku 
DO  26  L  =  1  *  LV 

S1=51+Z(K*L)*P(  i  *  K  )  *'J  (  J  *  L  ) 

2b  S2  =S2  +  P ( I  *  K ) *R ( 1 *R ) *0 ( J  * L / *Q ( J  * L ) 

A (  I  * J )  =  Si /S2 


*  ' 


•  ( 


*  'l 


. 


•  *  A 

*  = 

• 

• 

* 

.  .  x  -  )  C 

9 

• 

. 


V  u 

J  9 

<  J  •  \ 

♦  1 

*  *  ) 

f 

♦ 

J  # 

1  I  * 

• 

• 

* 

/' r  • 

.  1 

. 

• 

* 

• 

l  -  * 

- 1  *  1 

£-  — 


1 


2  5 
2  I 


4J 


^2 

4  1 
4  3 


4  1 

4U 


4  4 
4/ 


wtviih  (6  *  2  /  )  l  *  J  *  A  I  i  «  j  ) 

FORMAT  (  1  H  *  3  A  9  i  2  *  ^  A  *  I  2  *  '■  X  9  L  i  •  n  ) 

T  HF  ArrROX  I  vi  A T  I  Oin# 

W  R  I  T  F  (  A  9  v  j  ) 

FORMAT  (  1  H o  9  4 A  9  a f h a  (  K  )  V  (  L.  )  7  (  K  9  L  ) 

no  4i  i  =1  9  LV 

1)0  4  I  K=  1  9  KU 

Fxr=u,  0 
i)0  ^2  i  =  i  9  in  a 

DU  -S2  J  =  1  9  M  A 

p  a  Y  =  p  a  r  +a (  l 9 J ) (  J  *  A  ) *0 i j  9  l ) 

Ut.V  =  AHSirX'r-4(K9L)  ) 

w  R  I  T  t  (  6  »  4  4  )  aim  9  G  r  (  L  )  *  Zip’Ll  9  Fad  u  f  V 

FORMAT  ( 1  H  *  Dt  F  6  •  2  »  4  a  *  F  6  •  2  9  4X9  i  .  *  r  *  4  a  *  9*6  9  3X9  F I 

CONVFRSION  OF  ORTl  ML  TO  kOwI  hOI  IAL5. 

R  (  1  9  1  )  =  1 . 0  0 
R  (  2  9  2  )  =  R  l  I  9  I  ) 

R  1  4  9  I  )  =  -  A  X  (  1  ) 

PI  4  9  3 ) =  R l 7  9  2  ) 

H(  'GO=Ri/»l  )  -  AX  (  X  )  *P  (2*2) 

R (  4  9  1  )  =  — R (2*1) *Aa ( 2  ) -HA l 2  ) 

1)0  40  I  =4  9  IMA 
1  A  =  1  -  i 
JA= 1-2 

R  (  1  9  I  )  =  R  i  1  A  9  i  A  ) 

r  (  J,1A)=R(  I  A*JA  l-AAl  1  A  )  -*R  (  i  A  9  i  A  ) 

1)0  41  J  =  2  *  JA 
J  H  =  I-J 

R(  J  9  J  R  )  =  P  (  1A9JH  —  1  >  —  A  X  (  1  A  i  ■*  R  i  i  A  9  J-  1-hXi  IAi*P(  JA*J'  ) 
h l  ) = -A  X (  I  A  )  *R (  I  A  * ]  )  -hX (  I  A ) ( J a  *  ) 

u  (  1  9  1  )  =  ]  •  0  0 
0 (  7  9  2 ) =0 (  1  *  1  ) 

0(2*1  ) = -A  Y ( 1 ) 

Ui  4 9  3 ) =  0 ( 2  9  2 ) 
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